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FOREWORD 


This report documents the results of an effort to develop analytical methods for 
determining flow fields due to stratified and closely spaced jets exhausting into a cross- 
flow. In conjunction with this study, a wind tunnel test program was conducted to gen- 
erate data against which analytical results could be compared. 

Mr. M. F. Schwendemann directed the experimental phase of this investigation, 
which is documented in Northrop report NOR 73-98. 

The work was performed by the Northrop Corporation under NASA contract 
NASI -11524, under the technical cognizance of Mr. Richard J. Margason. 

This report has been assigned the number NOR 73-77 for internal control 
purposes. 
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ANALYSIS OF STRATIFIED AND CLOSELY SPACED JETS 
EXHAUSTING INTO A CROSSFLOW 

By 

H. Ziegler and P. T. Wooler 
Northrop Corporation, Hawthorne, California 

SUMMARY 


Procedures have been developed for determining the flow field about jets with 
velocity stratification exhausting into a crossflow. Jets with three different types of 
exit velocity stratification have been considered, namely 

• Jets with a relatively high velocity core 

• Jets with a relatively low velocity core 

• Jets originating from a vaned nozzle 

The procedure developed for a jet originating from a high velocity core nozzle 
is to construct an equivalent nozzle having the same mass flow and thrust but having 
a uniform exit velocity profile. Calculations of the jet centerline and induced surface 
static pressures have been shown to be in good agreement with test data for a high 
velocity core nozzle. 

The equivalent ideal nozzle has also been shown to be a good representation for 
jets with a relatively low velocity core and for jets originating from a vaned nozzle 
in evaluating jet-induced flow fields. 

For the singular case of a low velocity core nozzle, namely a nozzle with a dead 
air core, and for the vaned nozzle, an alternative procedure has been developed. The 
internal mixing which takes place in the jet core has been properly accunted for in the 
equations of motion governing the jet development. Calculations of jet centerlines and 
induced surface static pressures show good agreement with test data for these nozzles. 

A method for treating two-jet configurations, formulated in an earlier investi- 
gation, has been extended to include mutual interference effects between the two jets 
in addition to the jet blockage effects already considered. Comparisons are made 
between calculations and test data for a number of jet configurations. 
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INTRODUCTION 


A fundamental problem in the development of methods for predicting aerodynamic 
characteristics of lift-jet, vectored thrust and lift -fan V/STOL aircraft is that of for- 
mulating a mathematical model to estimate the effects of the propulsion system efflux 
interaction with a crossflow. During the transition flight phase, this efflux is directed 
at large angles to the freestream and has a significant influence on the aircraft aero- 
dynamics as well as on the stability and control requirements. Consequently, a con- 
siderable amount of research activity, both experimental and analytical, has been de- 
voted to the development of an understanding of this flow problem and also to the deve- 
lopment of methods to calculate the resulting interference flow fields. 

A number of analytical formulations of the problem of a single jet exhausting into 
a crossflow exist, and details of the different approaches may be found in reference 1. 
An approach to the problem of a single, normally exhausting jet, which appeared to 
offer possibilities of treating more complex flow configurations, may be found in 
reference 2. An entrainment model was developed from dimensional analysis and 
physical considerations. The force on the jet boundary, as a result of the pressure 
differential around the jet, was accounted for by a crossflow drag. The geometry of 
the jet cross section was represented by an ellipse and the continuity and momentum 
equations were solved to provide the jet path. The jet -induced velocity field was then 
determined by replacing the jet by a distribution of sinks and doublets. Using an image 
system for the flat plate and lifting surface theory for the finite wing, it was then 
possible to determine the jet-induced pressure distribution on these two types of sur- 
faces. 

Induced surface static pressure distributions around single jets exhausting nor- 
mally into a crossflow have been determined in references 3, 4, 5 and 6. Jet de- 
cay investigations for single jets exhausting at 90° into a crossflow were conducted in 
references 7 and 8. Jet centerlines (usually defined as the position of maximum total 
head in the jet) have been obtained for single jets exhausting at various angles into the 
crossflow (references 9, 10, and 11). 

The analytical model described above has been further extended in reference 12 
to treat jets exhausting into arbitrarily directed crossflows as well as multiple-jet 



configurations. In the case of multiple-jet configurations it was assumed that the 
leading jet (jet farthest upstream) develops independently of the downstream jets until 
intersection occurs. It was also assumed that downstream jets behave as single jets 
developing into a crossflow of reduced dynamic pressure. An arbitrary jet configuration 
could then be treated as a combination of discrete jets, with the induced velocity com- 
ponents due to each jet being additive at each control point. Data from the wind tunnel 
investigation of reference 13 were used to substantiate the assumptions made in the 
development of the analytical model. Empirical relationships postulated in the develop- 
ment of the model were established quantitatively. 

hi general, the exit velocity profile of a V/STOL fan or jet will not be uniform 
(see reference 14, for example). Variations in dynamic pressure decay may also 
exist. In the hover case, it has been shown (reference 15) that jets with different de- 
cay rates produce different induced aerodynamic forces, so that similar effects are 
expected with the jet exhausting into a crossflow. 

It has been the purpose of this study to develop methods for determining flow 
fields due to stratified and closely spaced jets exhausting into a crossflow. In con- 
junction with this study a wind tunnel test program has been conducted to generate 
data against which analytical results can be compared. 
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m ieff 
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q« 

q e 
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ratio of inner to outer diameter for an annular nozzle (see sketch 2) 

exit area of an equivalent ideal nozzle 

jet cross-sectional area 

core mixing parameter, ado / 1 

circumference of jet cross section 

crossflow drag coefficient of jet cross section 

pressure coefficient, (p-p«>)/qoo 

length of major axis in elliptical representation of jet cross section 

jet exit diameter 

d/do 

effective jet exit diameter, obtained from equivalent ideal nozzle 
considerations (see equation (11)) 

diameter of downstream jet (see sketch 4) 

entrainment crossflow per unit length of jet (see equation (1)) 

entrainment parameters (see equation (1)) 

nozzle thrust 

jet blockage factor, (d 0 -A)/do 

length of core region for an annular nozzle (see sketch 2) 

jet exit velocity parameter, inverse velocity ratio U j 0 /Uoo 

inverse velocity ratio for i^ segment of jet in mutual interference 
computations 

jet mass flow 

jet mass flow at exit 

Mach number 

static pressure 

crossflow static pressure 

jet dynamic pressure at exit of an ideal nozzle 

crossflow dynamic pressure, 

effective crossflow dynamic pressure for downstream jet when jets 
are aligned in freestream direction 

q e when leading jet does not exhaust normally into the crossflow 
(see equation 22) 
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SYMBOLS (Continued) 
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Qo 


effective crossflow dynamic pressure for downstream jet when jets 
are not aligned in the crossflow direction 

jet dynamic pressure, jet dynamic pressure at exit 

effective jet dynamic pressure at exit 

gas constant 

spacing between two jets 

stagnation temperature 

crossflow speed 

jet speed, jet speed at exit 

nondimens ionalized jet speed, Uj/Uj Q 

velocity ratio, [qoo/qj 0 ]' /j 

induced velocity vector for freestream modification (see sketch 5) 
Cartesian coordinate system 

distances in and normal to the crossflow direction normalized by do 

sideslip angle 

ratio of specific heats 

overlap of downstream jet by upstream jet 

jet deflection angle at exit (see sketch 8) 

density 

angle between local velocity vector and the normal to the crossflow 
vector (see sketch 1) 

value of 0 at jet exit 


Subscripts 

i value in i^ 1 segment of jet in mutual interference computations 

°i initial value for i tfl segment of jet in mutual interference computations 
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VELOCITY STRATIFICATION EFFECTS 


In general, the exit velocity profile of a V/STOL fan or jet engine will not be uni- 
form. Variations in dynamic pressure decay may also exist. The jet exit velocity 
profile and variations in dynamic pressure decay have an effect on the jet-induced ve- 
locity field in a crossflow and consequently influence the aircraft aerodynamics. 

In this report three different types of nozzle flows are considered, namely 

1. Annular nozzle with high velocity core 

2. Annular nozzle with low velocity core 

3. Vaned nozzle 

The annular nozzle with high velocity core is representative of the exhaust flow 
from a turbofan engine. The annular nozzle with low velocity core models the flow 
from a lift-fan engine. The vaned nozzle is sometimes considered for vectored thrust 
concepts. 

In earlier investigations, references 2 and 12, a mathematical model was devel- 
oped for the flow about a subsonic turbulent jet exhausting at an angle into a uniform 
crossflow. In these studies the jet was assumed to be deflected in the crossflow di- 
rection due to entrainment of crossflow fluid and also due to jet blockage. 

Consider a circular jet exhausting at a right angle into a uniform mainstream, 
as shown in figure 1. The entrainment of crossflow fluid was represented, in refer- 
ence 2, by the expression 


pu!d # 


E , d* cos 6 


E 2 (m U* - sin 0)C/d e 

+ 1 + E, cos 0/mU* 

3 J 


(1) 


The entrainment parameter E 2 is obtained from static results for the jet and 
may vary from nozzle to nozzle, whereas E-^ and Eg have been determined empirically 
and are independent of the nozzle configuration. 
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The equations of motion for the jet are, according to reference 2, 
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where primes denote differentiation with respect to Z*. 

From equations (1) through (4) it is observed that m is a parameter and that do 
is the length dimension for normalizing purposes. In the case of a standard conver- 
gent subsonic nozzle, m and d Q are clearly defined as the square root of the ratio of 
the jet exit to mainstream dynamic pressures and the nozzle exit diameter, respec- 
tively. One of the approaches to treating jets with stratified exit flow characteristics 
explored as part of this study is to employ the analytical model of reference 2, con- 
tingent upon an appropriate determination of m and d Q for the nozzle. 

In order that the analytical methods developed in this study may be used with 
confidence to calculate induced flow fields, pressure distributions and forces and 
moments on adjacent aerodynamic surfaces, it is desirable to make comparisons 
between calculations and test data for induced surface pressures on a simple geomet- 
ric shape due to stratified jets exhausting into a crossflow. 

To generate test data against which analytical results can be compared, a wind 
tunnel test program has been conducted and documented in detail in reference 16. 

A four-foot diameter circular plate containing the nozzles was aligned with the tunnel 
flow. The plate contained pressure taps to determine the surface static pressure 
distributions around the exhausting jets. Jet centerlines and decay characteristics 
were obtained from a total head rake. To obtain the three types of stratified exit 
flow characteristics considered in this study, the nozzles of figure 2 were utilized. 
The air supply for the core and annular regions of the dual concentric nozzle could be 
controlled independently to yield an annular jet with a high or a low velocity core. 
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Equivalent Ideal Nozzle 


The procedure which has been developed for determining an effective m and cIq 
for nozzles with stratified exit flow characteristics is to consider a nozzle of the same 
mass flow and thrust, but having a uniform exit velocity profile. 

The mass flow m for an ideal nozzle of area A is 


m p m (l + 2=L M f 


(5) 


where p and M are the static pressure and Mach number of the flow at the nozzle 
exit, respectively. 

The thrust F, assuming subsonic flow, is 

F = 2Aq 


( 6 ) 


where q is the nozzle exit dynamic pressure. Equation (5) may be written in terms of 
the dynamic pressure, so that 




(7) 


Eliminating A between equations (6) and (7) we deduce 



and then, substituting for q in equation (6) , we obtain 



( 8 ) 


0 ) 
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Thus, q and A for the ideal nozzle may be determined by substituting the nozzle mass 
flow and thrust in equations (8) and (9) . The parameter m is then determined from 
the expression 



( 10 ) 


q<x> being the dynamic pressure of the crossflow. The diameter d Q used for normalizing 
purposes is determined from equation (9) . It is 



Annular Nozzle with High Velocity Core 

By considering an ideal nozzle, of the same mass flow and thrust, the effective 
jet dynamic pressure and effective diameter for the high velocity core nozzle tested 
in the related wind tunnel investigation (reference 16) have been determined. The core 
of the nozzle had an area of 0. 8026 sq cm (0. 1244 sq in), the annular region had an 
area of 4.0948 sq cm (0. 6347 sq in). The dynamic pressures of the core and annular 
regions were 7.0931 N/cm 2 (1481.36 psf) and 2.7409 N/cm 2 (572.42 psf), respectively. 
Thus, assuming the jet total temperature to be ambient and the jet to exhaust into a 
standard atmosphere at sea level, we obtain 

qj e = 3. 5824 N/cm 2 (748. 18 psf) 
d e = 2.4518 cm (0.9653 in) 

The effective diameter is not too different from the outer diameter of the nozzle, 
which is 2. 54 cm (1.0 in), whereas the effective jet exit dynamic pressure is signifi- 
cantly different from both the core and annular values. The entrainment character- 
istics for the high velocity core nozzle, shown in figure 3, are not significantly differ- 
ent from those for the other stratified nozzles or for a convergent nozzle. 
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A sketch of the single- jet configuration tested in this investigation, showing the 
spanwise pressure tap stations at which comparisons between test data and calculations 
have been made, is shown below. 



SKETCH 1 

The theoretical prediction for the jet centerline is shown in figure 4. This may 
be compared with the test data for the high velocity core nozzle centerline shown in 
figure 5. Calculations of the surface pressure distribution compared with test data 
for the high velocity core nozzle, as well as for a clean nozzle at the same velocity 
ratio, are shown in figure 6. 

It should be pointed out that, while comparisons between analytical results and 
experimental data are shown for only a representative number of jet configurations 
(and generally restricted to two spanwise stations for each configuration), this should 
not be viewed as a limit on the range of data acquired in the experimental phase of 
this investigation. Complete and extensive documentation of the wind tunnel test 
program may be found in reference 16. 
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Annular Nozzle with Low Velocity Core 


The singular case of an annular nozzle with a low velocity core, namely a dead 
air core, has been investigated. The dual concentric nozzle of figure 2 was utilized 
with no air being supplied to the core region. A schematic representation of the 
velocity profiles at the exit and at the end of the core region is shown in sketch 2. 
Details on the actual velocity profiles and the decay of the jet dynamic pressure may 
be found in reference 16. 



t-T 


SKETCH 2 

After establishing the effective jet exit dynamic pressure for the nozzle, the 
tunnel (or freestream) dynamic pressure was adjusted to give the appropriate jet exit 
velocity ratio, and induced surface static pressure data were acquired. Figures 48 
through 50 of reference 16 show that the experimental data for induced surface static 
pressures around the annular jet with dead air core collapse on the data for a uniform 
jet of the same velocity ratio very well when the effective jet exit diameter of the 
annular nozzle is used in nondimensionalizing distances associated with those data, 
except in the wake region behind the jet. 
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Thus, the concept of utilizing an effective m and d Q , obtained from considering 
an equivalent ideal nozzle as discussed previously, in the analytical model of reference 
2 is seen to be a valid approach to treating velocity stratification effects for the annu- 
lar nozzle with low velocity core in computing induced surface static pressures around 
the jet in those regions where the model of reference 2 is expected to yield meaning- 
ful results. The singularity representation of the jet used to compute the jet-induced 
pressure distribution in reference 2 assumes potential flow in the areas external to 
the jet and constitutes a good representation only in those regions where the flow 
outside the boundary layer is potential. 

As an improvement to the equivalent ideal nozzle approach, the analysis of 
reference 2 may be modified for the annular jet with dead air core so that the internal 
mixing which takes place in the jet core may be properly accounted for. 

Consider the annular jet of sketch 2. Let the diameter of the core, which is a 
region of no efflux (qj 0 = 0), be ado, and the outer diameter of the annulus be do, and 
let the length of the core region be i. The diameter of the core at distance Z from 
the nozzle is ad Q -bZ, where b = ad 0 /i . The equations of motion for a single jet ex- 
hausting into a crossflow are given by equations (2) through (4). 

In the original analysis of uniform turbulent jets originating from circular noz- 
zles, a development region was assumed in which the jet deformed from a circular 
cross section into an elliptical one. It is assumed here that this development region 
is not changed by the core region and that, in turn, the core region is independent of 
crossflow velocity. 

Thus in the core region, 0 ^ Z* g a/b, we obtain 

di = T d * 2 [ 1 ' TS ' d * 5 ( a ' ” 2ab z * + b2z * 2 )] (12) 

for the development region, Z*^ 0. 3m, and 

^ = t£ ! [i - 3^ (a ! - 2ab Z* + bM] < 13 > 


for the developed region, Z* > 0. 3m. Outside of the core region, Z* >a/b, the ex- 
pressions for Aj/d 0 2 are as above with a = b = 0. 
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The expressions for the jet circumference, C, are not changed by the inclusion 
of the core term, so that 
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Equations (2) through (4) then become 
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where the parameter m for this nozzle is now defined to be 
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Equations (16) and (17) are valid for both the development region and the devel- 
oped region of the jet. Equation (18) is to be used for the development region of the 
jet and equation (19) is valid for the fully developed region. Outside the core region, 

Z* >a/b, equations (18) and (19) are modified by setting a = b = 0. 

Equations (1), (12) - (19), together with the jet exit boundary conditions, are the 
equations governing the jet development. They may be integrated, following reference 
12, to obtain UJ, d* and X* as functions of Z*. 

The annular jet which was tested in this investigation had an outer diameter of 
2. 54 cm ( 1. 0 in) and a core diameter of 1.02 cm (0.40 in). The jet exit dynamic 
pressure was 7.09 N/cm 2 (1481 psf) for the annular region and the core region ex- 
tended over a length of 10. 16 cm (4.0 in). The parameters a, b are then 0.40 and 
0.10, respectively. 

The computed jet centerline for the annular nozzle at a velocity ratio of 0. 125 
is shown in figure 4. The corresponding test data are shown in figure 5. From these 
results it is observed that the annular jet does not penetrate into the crossflow as far 
as a high velocity core nozzle (or clean nozzle) at the same velocity ratio. 

Surface static pressure calculations are shown in figure 7, with the test data 
also shown for comparison. Good correlation between theory and test data is observed. 
Comparison with figure 6 indicates that the annular nozzle induces surface pressures 
of slightly smaller magnitude than the high velocity core nozzle (or clean nozzle). 


Vaned Nozzle 

The vaned nozzle which has been considered in the experimental phase of this 
investigation is shown in figure 2. The induced surface static pressure data plotted 
in figures 48 and 49 of reference 16 (again utilizing for nondimensionalizing purposes 
the effective jet exit diameter for the vaned nozzle, obtained from equivalent ideal 
nozzle considerations) indicate that the induced pressure distributions around a vaned 
jet may be predicted quite accurately by using the equivalent ideal nozzle approach. 

The presence of the vanes in the nozzle reduces the exit area and one might 
expect a change in the mixing characteristics for this nozzle. However, if the mass 
flow m is plotted against distance, normalized by the effective diameter as defined in 
equation (11), the entrainment characteristics of all the nozzles are similar (figure 3). 
It may be deduced, therefore, that the vanes only affect the internal mixing of the jet. 
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so that the modified approach developed for the annular nozzle with dead air core, to 
account for the internal mixing taking place in the jet core, may serve as an improve- 
ment over the equivalent ideal nozzle approach for the vaned jet as well. The area of 
the vanes at the nozzle exit is replaced by a circular core of equal area and the extent 
of the core region is again determined from static test results (reference 16). 

The nozzle which was tested in this study had a vaned area of 1. 719 sq cm 
(0.266 sq in), yielding a core diameter of 1. 48 cm (0.58 in). The extent of the core 
region was determined to be 14. 73 cm ( 5. 8 in). The parameters a and b of equation 
(12) are, therefore, 0.58 and 0. 10, respectively. The parameter m for this nozzle 
was defined as 


in 



mainstream dynamic pressure 
jet exit maximum dynamic pressure 


Computed jet centerlines for the vaned nozzle, at velocity ratios of 0. 125 and 
0.250, are shown in figure 4. The corresponding test data are shown in figure 5. 
Figure 5 includes test data for the vanes perpendicular to the crossflow and for the 
vanes aligned with the crossflow. The calculations of figure 4 do not, of course, 
account for this difference in the alignment of the vanes. The vaned nozzle jet is ob- 
served to penetrate the mainstream less than either the annular nozzle jet or the high 
velocity core jet at the same velocity ratio. 

Figure 8 shows computed surface pressure distributions for the vaned nozzle, 
with the vanes perpendicular to the crossflow, at a velocity ratio of 0. 125. Test data 
are included for comparison. The correlation is observed to be quite good. The pres- 
sure data for this nozzle, both calculated and test, are observed to be of slightly 
lower magnitude than those for the annular nozzle (figure 7), at the same velocity ratio. 
Computed and experimental surface pressure distributions for the vaned nozzle, 
with the vanes aligned with the crossflow, are shown in figure 9. The velocity ratio 
is again 0. 125 and the computed pressure distributions are, of course, the same 
as those shown in figure 8. No significant effect due to the changing of the orientation 
of the vanes is observed from the test data. 

Figure 10 shows pressure distributions around the vaned nozzle jet, with the 
vanes perpendicular to the crossflow, at a velocity ratio of 0. 250. The test data from 
two runs using different values of jet exit dynamic pressure and freestream dynamic 
pressure to achieve the velocity ratio are included. No significant differences between 
these two sets of data are observed. 
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MUTUAL INTERFERENCE EFFECTS 

The single-jet model of reference 2 was extended in reference 12 to treat jets 
exhausting into arbitrarily directed freestreams, as well as multiple-jet configurations. 

Multiple-jet configurations were treated as combinations of discrete jets, with 
leading jets assumed to develop independently and downstream jets assumed to exhaust 
into a crossflow of reduced dynamic pressure. The test data of reference 13 were 
used to derive an expression for the effective dynamic pressure which the downstream 
jet "sees" as a result of the blockage of the crossflow by the upstream jet. This ex- 
pression was based on data for two-jet configurations at zero sideslip (jet exits aligned 
in the freestream direction), with both jets exhausting normally into the freestream. 

The test data were utilized to verify that the upstream jet develops independently 
of the downstream jet for the zero sideslip condition, even for the closely spaced two- 
jet configuration tested in reference 13. For the closely spaced configuration at zero 
sideslip, the downstream jet "sees" a low crossflow dynamic pressure and, conse- 
quently, does not exert a strong influence on the induced flow field. The assumption 
that the upstream jet develops independently of the downstream jet is therefore justi- 
fied, despite the close jet spacing. Since the expression for the effective crossflow 
dynamic pressure for the downstream jet was based on data for the zero sideslip con- 
figurations, it accounts for all the interference effects between the two jets. The good 
agreement between calculated induced pressure distributions and the test data for the 
zero sideslip configurations, exhibited in the comparisons of reference 12, supports 
this conclusion. 

For the closely spaced configuration at sideslip, noticeable differences between 
theory and test data were evident. With the jets no longer aligned in the freestream 
direction, the downstream jet now has a stronger influence on the induced flow field 
since there is less blockage of the crossflow by the upstream jet. This stronger in- 
fluence, together with the close jet spacing, makes the assumption that the upstream 
jet develops independently of the downstream jet no longer representative of the 
physical situation. It was felt that further mutual interference effects between the two 
jets had to be included to improve correlation between theory and test data. 
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The extension of the basic two-jet computation to include mutual interference 
effects between the two jets, in addition to the blockage effect discussed previously, 
is presented here. Comparisons between theory and test data of reference 13 for the 
closely spaced two-jet configuration are shown. 

The modification of the expression for the crossflow dynamic pressure which the 
downstream jet "sees" in a two-jet configuration, with both jets exhausting at an angle 
other than 90° into the crossflow, is discussed. Calculations of jet centerlines and in- 
duced surface static pressures are compared with test data of reference 16 for three 
different spacings between two inclined jets . 


Two- Jet Analytical Model 

The details of the computational procedure applying the basic single-jet model 
to the calculation of the interaction flow field due to two exhausting jets are given in 
reference 12. 

A two-jet configuration was treated as a combination of discrete jets. The equa- 
tions of motion for each of the exhausting jets were integrated, utilizing the appropriate 
initial conditions for each jet, Z* = 0. , X* = 0. , Uj* = 1 . , d* = 1 . , and, using the coor- 
dinate system of figure 1, dX*/dZ* = [(1. -cos 2 0 o )/cos 2 0 o ]’ /2 , as well as the correspon- 
ding jet exit velocity ratio, to yield the mean jet speed U?, the major diameter of the 
ellipse representing the jet cross section d*, and the displacement of the jet center- 
line in the freestream direction X*, all in nondimensionalized form, as functions of 
Z*, the nondimensionalized penetration of the jet centerline into the crossflow. 

The upstream jet was assumed to develop independently of the downstream jet 
and the downstream jet was assumed to exhaust into a crossflow of reduced dynamic 
pressure, which it "sees" as the result of blockage by the upstream jet. Thus the in- 
fluence of the upstream jet on the downstream jet was introduced into the computations 
as a reduced freestream velocity, Ug/Uoo = [ q e / ] ^ 2 , in the equations governing the 
development of the downstream jet ( equations (2) - (4) ) . 

The extent of overlap between the two jets was the principal parameter in deter- 
mining the degree of blockage experienced by the downstream jet. The computational 
details of establishing this degree of influence of the upstream jet on the downstream 
jet for each jet element, as the integration of the equations of motion is being carried 
out, are given in reference 12. 
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Test data of reference 13 were used to obtain an empirical relationship for the 
dynamic pressure q e which the downstream jet "sees" as the result of the crossflow 
blockage by the upstream jet, in terms of the crossflow dynamic pressure, qoo, and 
the spacing between the two jets, s (see sketch below). 


Jet# 1 Jet #2 



SKETCH 3 

This expression is given in reference 12 as 



s/d 0 - 1 
s/dp + .75 


{s/dp > 1} 


( 20 ) 


and is used as a limiting value, when computation of the overlap between the two jets 
shows that the downstream jet element is completely in the zone of influence of the up- 
stream jet element (as, for example, in the case of two jets aligned in the crossflow 
direction) . When the two jets are not aligned, an effective crossflow dynamic pressure, 
q^, which is a weighted mean of q e given above and q^,, is utilized. The weighting of 
the dynamic pressure is determined from the degree of overlap between the upstream 
and the downstream jet elements discussed previously and shown in schematic form 
below. 



SKETCH 4 
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Thus, 


^ * (d,-A)V^ pi) 

d 2 

As the equations of motion for the two exhausting jets are being integrated, the 
distance between the two jet centerlines is continually checked. When intersection of 
the two jets is indicated, initial conditions for the merged jet which results are deter- 
mined from continuity and momentum considerations, as detailed in reference 12. 

These initial conditions are then employed in integrating the set of differential equa- 
tions for u| , d*, and X*. 

The velocity field induced by a two-jet configuration can now be determined by 
replacing each jet (including the jet resulting from the coalescence of the two exhaust- 
ing jets) by its representative singularity distribution of sinks and doublets . 

The expression for the dynamic pressure to be utilized in the downstream jet 
computations (equation (20)) was based on data for two-jet configurations at zero side- 
slip, with both jets exhausting normally into the crossflow, and contains only the jet 
exit spacing as a parameter. 

Since equation (20) was based on data for zero sideslip configurations, it account- 
ed for all the interference effects between the two jets, even for the closely spaced 
configuration, s = 2. 5d 0 , of reference 13. Comparisons of test data with theoretical 
predictions in reference 12 showed good correlation for the zero sideslip configurations 
(see, for example, figures 37 and 40 of reference 12) . 

Noticeable differences between theory and test data were discernible for the 
non-zero sideslip configurations, particularly for the close jet spacing (see, for ex- 
ample, figure 41 in reference 12). For these computations q e as given by equation (20) 
is weighted with q^ according to equation (21). The downstream jet now has a stronger 
influence on the induced flow field, and it was felt that, although the relationship of 
equation (21) accounts for the effect of blockage of the upstream jet on the downstream 
jet, further mutual interference effects between the two jets had to be included in the 
computations to improve the correlation between theory and test data. 

For the purpose of including mutual jet interference effects, an iterative pro- 
cedure involving modification of the crossflow into which the two jets exhaust, has 
been chosen. The scheme treats each of the exhausting jets in a number of segments 
as shown in sketch 5. 
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The freestream velocity vector for each segment of one exhausting jet is per- 
turbed by the induced velocity vector (u, , u 2 , etc) due to the other exhausting jet. 

The perturbed freestream velocity vector is assumed constant over the extent of the 
segment and is evaluated at the point of origin of each segment. Each segment is 
treated as a discrete jet, with proper initial conditions and the appropriate freestream 
velocity vector. 

The first computation sets u, , u 2 , etc, equal to zero to establish the first approx- 
imation for the centerlines of the two exhausting jets and the coalesced jet, if inter- 
section between the two jets occurs. The two exhausting jets and the coalesced jet 
are then replaced by their representative singularity distributions and the induced 
velocities u, , u 2 , etc, are then computed for each segment of the two exhausting 
jets. 

For the first iteration, the initial conditions for segment I now become 
d? = l., Uj", = Uj,/Uj 0 , = 1. » m, = Uj 0l /Uoo 

where m, is the inverse velocity ratio of the exhausting jet. 

The direction cosines of the modified freestream velocity vector, U® +u,, are 
determined and a local coordinate system is established which is aligned with the free- 
stream velocity vector and the jet exhaust vector (see reference 12, page 25) . The 
initial value for dX*/dZ* is determined in this local coordinate system. 
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The effective inverse velocity ratio for segment I is 


m i eff = 


Ujo, / P» 

Uoo yUec+U, 



Uqo 

Uoo +U, 


) 


The equations for U*, d*, X* may now be integrated over the extent of segment I. 

The last point of segment I then becomes the origin of the next segment (or next 
discrete jet) with a diameter d 02 = d* d 0l , where d* is the last computed value of the 
nondimensionalized jet diameter in segment I 

Other initial conditions for segment n are 

d 2 =l., Uj 2 = Uj 2 /Uj 02 = 1. , m 2 =Uj 02 /Uoo =Uj,m, 


where Uj^ is the last computed value of the nondimensionalized mean jet speed in 
segment I. 

At this point the direction cosines of TTqo + u 2 are determined and a new local, 
jet-oriented coordinate system is established. The initial value for dX*/dZ* is deter- 
mined in this coordinate system from the known direction of the jet centerline at the 
end of segment I. The effective inverse velocity ratio for segment II is 


m 2eff = 



Uji m , 


Uoo 

Uoo+U 2 


) 


The equations for Uj , d*, X* may now be integrated over the extent of segment II. 

The computations described above continue for each of the two exhausting jets 
until integration over the extent of each jet, up to the point of intersection with the 
other jet, has been accomplished. The procedure of establishing intial conditions 
for the jet resulting from coalescence of the two exhausting jets remains unchanged 
from that detailed in reference 12 and outlined previously. New values for u,, u 2 , etc, 
may now be obtained and the entire computational scheme may be repeated. 

In determining the induced velocity vectors uj for the segments of each exhaust- 
ing jet, only the effect of the other jet is to be considered. In the representation of 
sketch 6, the singularity distributions of two exhausting jets, ©and (5), and of the 
coalesced jet© are indicated. The singularity distribution©) is a continuation of © 
and would result if no intersection occurred between Jet#l and Jet #2. 


20 



Jet#l 


Jet #2 



SKETCH 6 

It is desired to evaluate the flow field induced at point P by Jet #2 alone. This 
is accomplished by summing the contributions to the induced velocity components at P 
from the segments constituting (5) and adding the contributions due to the coalesced 
jet ( 5 ) But the coalesced jet is established from continuity and momentum consider- 
ations involving both Jet#l and Jet #2 and thus, some influence of Jet#l at point P is 
now included. The contributions to the induced flow field at P due to © are now sub- 
tracted to account for this, and the induced velocity components V x , Vy, V z are ob- 
tained. 

As discussed previously, the application of equation (20) in the two-jet computa- 
tions appeared to adequately account for all mutual interference effects for the zero 
sideslip configurations, since equation (20) was derived from zero-sidelip data. With 
increasing sideslip, equation (21) shows that q^ approaches q^ as A (or the degree of 
overlap) approaches zero. For a spanwise configuration, the only interference effects 
between the two jets would be those accounted for by the modification of the crossflow 
by the iterative procedure. 

Equation (21) then suggests that, after the induced velocity components V x , Vy, 

V z have been evaluated, the degree of overlap, as represented by the term (d 2 - A)/d 2 , 
be considered before the freestream velocity vector is modified. 
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As an approximation to the term (d 2 - A)/d 2 in equation (21), a factor K is eval- 
uated, using properties in the plane of the jet exits, as shown in sketch 7, such that 
K= (d 0 -A)/do. 



SKETCH 7 

The induced velocity vector is then obtained from the induced velocity 
components 

Si=K(V x t + V y t + V z fc) 

Thus, no further mutual interference effects are included for the zero sideslip 
configurations, and full effect of the modification of the crossflow is included for a 
configuration where equation (21) does not provide for any blockage effect. 


Comparisons of Two-Jet Calculations With Test Data 

Computations have been carried out primarily for the closely spaced two-jet 
configuration for which induced pressure distributions in the plane of the jet exits were 
obtained in reference 13. A schematic of this configuration is shown in sketch 8. 

Figure 11 shows a comparison between theory and test data for induced pressure 
variation with X/d 0 at Y/d Q = 1. 5 and 3. Since mutual interference effects due to the 
modification of the crossflow are not included for a zero sideslip configuration, as 
discussed previously, the computed pressure distributions should agree with those 
presented in figure 40 of reference 12. Comparison with the theoretical pressure dis- 
tributions of figure 40, reference 12, does confirm that numerical differences incurred 
by breaking the two exhausting jets into segments and treating the segments as discrete 
jets are negligible. This, of course, applies only to the computation of the initial 
approximation when Q,, U 2 , etc, are equal to zero. 

Figure 12 shows the same comparison for a sideslip angle /3 = 20°. Pressure 
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. . . . Pressure 
Taps 



SKETCH 8 

distributions corresponding to the first approximation and the final iteration of the 
computational scheme described previously are shown. Some improvement in corre- 
lation between theory and test data is discernible due to the incorporation of these 
further mutual interference effects. 

For this configuration, the projections of the centerlines in the X-Y and the 
X-Z plane are shown in figure 13. In contrast to the results for 13 = 0 (reference 12), 
the calculated centerline for the merged jet indicates greater penetration than observed 
in the wind tunnel test. Also, there is a significant difference between calculation and 
test data in the projection of the centerlines in the X-Y plane. The experimentally de- 
termined centerlines (positions of maximum total head) exhibit an unexpected lateral 
deflection, which may be due to the partial blockage of the downstream jet resulting in 
the jet momentum decaying at a decreased rate on the side of the jet which is blocked 
from the crossflow by the leading jet. 

A two-jet configuration with a spanwise spacing of 2. 5 diameters is shown in 
sketch 9. Theoretical and experimental pressure variations are shown for three 
stations of constant X/d 0 in figure 14. Again, the first approximation and the final 
iteration for each computed pressure distribution are shown. The full effect of the 
mutually induced velocities on the crossflow is included in the iterative procedure for 
this configuration, resulting in noticeable improvement in the correlation between 
theory and test data. 
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SKETCH 9 

The projections of the jet centerlines for this configuration, as well as for one 
with a spacing of 7. 5 diameters between the jet exits, are shown in figure 15. For the 
wider spacing, the jets are attracted to each other and each is deflected by the cross- 
flow to a greater extent than a single jet of the same velocity ratio. These features 
are predicted quite well by theory. For the close spacing, the analytical model pre- 
dicts that each jet will be deflected to a greater extent (up to the point of intersection) 
than the jets with the wide spacing. Following intersection, this trend is reversed 
and the computed centerline of the merged jet exhibits greater penetration than the 
individual jet centerlines of the widely spaced configuration. As was the case for the 
closely spaced configuration at sideslip j8= 20°, the positions of maxim um total head 
show an unexpected lateral deflection, in contrast to the calculation indicating inter- 
section after a penetration of about 6 diameters. 

In computing the pressure distributions of figures 12 and 14, two iterations on 
mutual interference effects were employed after the initial approximation to each 
pressure distribution had been established. Indications are that two iterations are 
sufficient. Experience with the test cases has shown that the first iteration (estab- 
lishing the induced velocities by which the crossflow is modified as other than zero) 
is primarily responsible for effecting the changes in the computed pressure distribu- 
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tions, and subsequent iterations produce little change. This is illustrated in figure 16 
where the first and second iterations for the pressure distributions shown in figures 
12 and 14 are presented. 


Multiple Inclined Jets 

The blockage effect of the upstream jet on the downstream jet, as given by 
equation (20), is seen to be a function of the spacing between the two exhausting jets 
only. The expression was based on experimental data of reference 13, where all 
multiple-jet configurations tested consisted of jets exhausting normally into the cross- 
flow. 

It was felt that for an inclined jet exhausting into the crossflow, equation (20) 
represents a reduction in the crossflow velocity normal to the jet. There is also a 
component of the crossflow velocity tangential to the jet, as shown below. 



SKETCH 10 

Thus, the magnitude of the effective dynamic pressure is 

= si"’0o] /! (22) 


where is given by equation (20). 

Figures 17 through 23 show comparisons of theoretical predictions and test data 
from reference 16, for three different spacings of two inclined jets. Equation (22) was 
utilized to account for the effect of blockage on the downstream jet and the iterative 
procedure described previously was employed for including further mutual interference 
effects between the two jets in the closely spaced configuration with sideslip. 
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The jet centerlines are predicted quite well. After intersection, the merged 
jet shows greater penetration into the crossflow than is indicated by theory, as was 
the case for jets exhausting normally into the crossflow (reference 12). The induced 
pressure distributions (figures 20 through 23) exhibit features similar to those dis- 
cernible in surface static pressure distributions due to two-jet configurations exhaust- 
ing normally into the crossflow. For the close spacing, the pressure distribution 
resembles that induced by a single jet. As the spacing between the two jets increases, 
the downstream jet is seen to have its own discrete effect on the pressure distribution. 
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CONCLUSIONS 


The concept of utilizing an effective jet exit velocity and diameter, obtained by 
considering a nozzle of the same mass flow and thrust but having a uniform exit veloc- 
ity profile (equivalent ideal nozzle), in conjunction with an analytical model for a uni- 
form jet exhausting into a crossflow has been shown to be a valid approach for deter- 
mining flow fields due to jets with exit velocity stratification. For the singular case 
of a low velocity core nozzle, namely a dual concentric nozzle with a dead air core, 
the analysis for a jet with uniform exit velocity profile has been modified to take into 
account the internal mixing in the jet core, as an improvement over the equivalent 
ideal nozzle approach. This modified analysis serves as an improvement over the 
equivalent ideal nozzle approach for the vaned jet as well. 

The investigation has shown that induced surface static pressure distributions 
due to stratified jets exhausting into a crossflow are not appreciably affected by the 
details of the exit velocity stratification. This indicates that small scale testing may 
"be accomplished with uniform exit velocity profile nozzles, without having to take re- 
course to reproducing, in detail, the stratified exit flow characteristics of lift/propul- 
sive systems of V/STOL configurations. 

Calculations of jet centerlines for the three types of nozzles investigated show 
that, for a given velocity ratio, the jet originating from a high velocity core nozzle 
penetrates the crossflow to the greatest extent, and the jet originating from the vaned 
nozzle exhibits the least penetration. These trends may be observed in the test data 
generated as part of this study. 

Inclusion of mutual interference effects between jets in a two-jet configuration 
has improved the correlation between theory and test data for two closely spaced jets 
exhausting normally into the crossflow. For two jets exhausting at an angle of 60 into 
the crossflow, calculations of centerlines and induced surface static pressures are in 
good agreement with test data. 
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FIGURE 8. INDUCED PRESSURE VARIATION FOR A VANED JET 
(Vanes Perpendicular to Freestream, Uoo/Uj o =0. 125) 
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FIGURE 11. INDUCED PRESSURE VARIATION FORA TWO-JET CONFIGURATION 
AT ZERO SIDESLIP (Spacing = 2. 5d 0 , <5j=90°, Uoo/U jo = 0. 125) 
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(b) Station Y/d 0 = 3 


FIGURE 12. (Concluded) 
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(a) Station X/d 0 = 0 


FIGURE 14. INDUCED PRESSURE VARIATION 

FOR A SPANWISE TWO-JET CONFIGURATION 
(Spacing = 2. 5 d 0 , <5j=90°, U<»/Uj 0 =0. 125) 
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FIGURE 21. INDUCED PRESSURE VARIATION FOR TWO INCLINED JETS 

( Spacing = 5 d 0 , <5i = 60°, U«/Ui o =0. 125) 
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FIGURE 22. INDUCED PRESSURE VARIATION FOR TWO INCLINED JETS 

(Spacing =7. 5 d 0 , <5j = 60°, U<»/Uj 0 = 0. 125) 


s 




Distance, X/d Q 

(b) Station Y/d 0 = 3 


FIGURE 23. INDUCED PRESSURE VARIATION FOR TWO INCLINED JETS 
AT SIDESLIP & =20° ( Spacing = 2. 5 d G , dj = 60°, U»/Uj o = 0. 125) 






APPENDIX A 


APPLICATION OF COMPUTER PROGRAM 
Sample Problem 

The use of the computer program in determining jet deflections and the jet- 
induced flow field is demonstrated for a general two-jet configuration exhausting into 
a crossflow. 

The two-jet configuration used in this sample problem is not one of the configu- 
rations considered in the experimental phase of this study. Rather, it is a composite 
configuration designed to exercise the computer program in its most general mode. 
This configuration is shown in figure Al. Nozzle #1 is the annular nozzle with dead 
air core descibed on pages 10-13. Nozzle #2 has the exit flow characteristics of the 
vaned nozzle described on pages 13 and 14, but has an exit diameter twice as large. 
Differing jet exit velocity ratios have been used for the two jets in the configuration. 

Input Data for Sample Problem 

The input data cards required for the sample problem are tabulated in figure A2 
and are described below. 

Card 1 lists five control indices. The first one, MULT=2, indicates that a two- 
jet configuration is being treated. The second one, IGE0M=4, specifies that pressure 
coefficients, as well as induced velocity components, are to be evaluated at all the 
control points provided as part of the input. By setting IPUNCH=0, no punched output 
is generated. The next two control indices deal with the calculation of mutual inter- 
ference effects, as discussed on pages 16-20. NPS=10 specifies that there will be 
10 integration intervals in each segment of each exhausting jet for which a modified 
freestream velocity vector will be computed. After establishing an initial approxima- 
tion using the unperturbed crossflow velocity vector, two iterations, utilizing cross- 
flow velocity vectors modified by mutually induced velocities for each segment, are 
specified by ITER=2. 

Card 2 specifies the angle of attack, a = 0, and angle of sideslip, (3 =20°, in the 
coordinate system of figure Al . 
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Card 3 controls the number of intervals and the interval size in the numerical 
integration of the equations of motion for the jet path. The integration is carried out 
in a local, jet-oriented coordinate system defined by the freestream velocity vector 
and the initial jet exhaust vector (for details, see reference 12). The numerical inte- 
gration routine in the program will optimize the actual integration step size utilized. 
Data for the jet centerline will be printed out at the specified intervals. For the 
sample problem, 80 intervals and an interval size of 0. 5 jet exit diameters are chosen. 

Cards 4, 5 and 6 describe the upstream jet. The jet location, in the coordinate 
system of figure Al, isX = -3. 75, Y = 0.0, Z=0.0. The jet exhaust angles (p and i// 
are 180° and 0, respectively. The jet exit diameter, d o = 1.0, and the jet exit velocity 
ratio, Uoo/Uj 0 =0. 125, are given. The parameters on card 6 serve to describe the 
exit flow characteristics of the nozzle. The ratio of effective core diameter to jet exit 
diameter for the annular nozzle with dead air core is 0. 40 and the jet mixing parameter 
is 0.10 (see discussion, page 13). 

Cards 7, 8 and 9 describe the downstream jet, which is located at X=-l, 25, 

Y = 0.0, Z - 0 . 0 , and again exhausts normally into the crossflow. The jet exit diameter 
is 2. 0 and the jet exit velocity ratio is 0. 250. The stratified exit flow characteristics 
of the vaned nozzle are specified by the parameters 0. 58 and 0. 10 (see page 14). 

Card 10 lists the parameter controlling the initial cross section of the jet re- 
sulting from the coalescence of the two exhausting jets. An ellipse with a minor to 
major axis ratio of 0. 5 is specified. (See Appendix B for guidelines on this para- 
meter) . 

Card 11 lists the number of spanwise stations, NS=3, and the number of control 
points at each station, NC=4, where the induced flow properties are to be evaluated. 

Cards 12-17 list the coordinates of the control points. All lie in the plane of the 
jet exits. The coordinates for each control point appear in the order X, Y, Z. The 
total number of control points is NCxNS. The listing is continuous, i. e. , no new rec- 
ord is required for the start of each spanwise station. 

Output for Sample Problem 

For the option specified on card 1, only printed output is obtained. Figure A3 (a) 
shows the first page of printed output. The jet configuration being treated is identified 
both by appropriate heading as well as other pertinent input data. Input controlling 
the numerical integration procedure is also displayed. 
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Figure A3(b) shows the jet centerline computations for the initial approximation 
(i. e. , the crossflow velocity vector is unperturbed). The coordinates of the jet cen- 
terline, the nondimensionalized mean jet speed, Uj/Uj Q , and the nondimensionalized 
major diameter of the ellipse representing the cross section of the jet, d/do» are 
printed out for each exhausting jet up to the point of coalescence. The point of coales- 
cence of the two jets is identified and initial conditions for the resulting jet are printed 
out. Jet centerline information for this jet, resulting from the intersection of the two 
exhausting jets, is then displayed. Jet centerline data are printed out at each inte- 
gration interval specified on card 3 of the input data, since for a normally exhausting 
jet at zero angle of attack the local jet-oriented coordinate system in which integra- 
tion is carried out and the fixed input/ output coordinate system coincide. The output 
in figure A3(b) displays only a portion of the jet centerline data generated for the 
sample problem. Computations would, of course, extend to Z = -40.0, which repre- 
sents integration of the equations of motion over the range |Z| = 80x0. 5x1.0 = 40.0, 
as specified on card 3 of the input data. 

Figure A3(d) shows the printout for the jet-induced velocity components and 
pressure coefficients at the control points specified. The coordinates of the control 
points are identified. The pressure coefficients and the induced velocity components 
U, V, W, nondimensionalized by Uoo, are given. 

The printout of figures A3(b) and A3(d) would then be repeated for the number of 
iterations specified on card 1. Figures A3(c) and A3(e) show the computations for 
the second, or final, iteration for the sample problem. Note that now, with the cross- 
flow velocity vector modified by mutually induced velocities for each s egment, each 
jet-oriented coordinate system is no longer aligned with the fixed input/output coordi- 
nate system. Thus, printout no longer occurs at the regular intervals of 0. 5, 1. 0, 

1. 5, etc. The program does adjust the local-coordinate integration interval to main- 
tain consistently spaced print on jet centerline data. 

Applicability and Limitations 

The program may be utilized to evaluate the induced flow field due to one or two 
jets exhausting into a crossflow. Jet exit velocity stratification effects may be treated 
by employing the velocity ratio and jet exit diameter for the equivalent ideal nozzle, 
or, alternatively, by accounting for the internal mixing through the introduction of 
the parameters describing jet exit flow characteristics. 

For a single-jet configuration, the initial jet exhaust direction, specified by <p 
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and if/ , and the freestream direction, specified by a and /3 , may be arbitrary. For 
a two-jet configuration, the jet exits must both lie in the same XY plane and the jet 
exhaust planes, defined by the freestream vector and the initial jet exhaust vectors, 
must be parallel. 

Comparisons between computations and experimental data have been made for 
velocity ratios 0. 10^Uoo/Uj o ^0. 30, and the program may be considered most appli- 
cable in this range. 

The choice of variables governing the numerical integration for the jet path is 
related to the velocity ratio of the problem being considered. For Uoo/Uj 0 ^ 0. 125, 
integration in the direction normal to the freestream over an extent of at least 30 
jet exit diameters has been found desirable. As U<»/Uj 0 increases, this may be de- 
creased, as the jet penetrates less at the higher velocity ratios. For the above range 
of velocity ratios an integration interval size of ^0.5 jet exit diameters has been found 
satisfactory. 

Control points at which jet-induced properties are to be evaluated may not lie 
within the jet itself, as the theory is not valid in this region. Generally, control points 
positioned less than 2 jet exit diameters from the center of a jet exit should be avoided* 
to minimize distortion in the computed velocity distributions. 

The jet -induced velocity field may be employed to explore the interaction between 
exhausting jets and adjacent supporting structures. Loading on adjacent lifting sur- 
faces may be evaluated by lifting surface theory, and other techniques, such as the 
transformation method, may be utilized for this purpose on fuselages. Details on the 
application of these methods may be found in reference 12, where the treatment of 
more complex multiple-jet configurations is also discussed. 
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(b) Jet Centerlines (Initial Approximation) 


FIGURE A3. (Continued) 


60 



*** ITERATION NUMBER ?t FINAL J TFp AT ION *<** 


** CFNTFRLTMFS OF JETS 1 AMD 2 ** 
AND rOALF^CEO JET 


X COORD 

YCOOPO 

ZCOOPO 

DJ 

ni a 

xcnoRn 

YCOORD 

7COORD 

U J 

D I A 

-3.75 

0*00 

0.00 

1 .000 

1 .00 

-1.25 

o.oo 

0.00 1 

.000 

1.00 

-3.75 

-.00 

-.50 

.945 

1.12 

-1.25 

-.00 

-.50 

.953 

1.09 

-3.74 

-.01 

-1.00 

.892 

1 .29 

-1.23 

-.01 

-1.00 

.908 

1.20 

-3.72 

-.01 

-1 .50 

.8^9 

1 .54 

-1.21 

-.02 

-1.50 

.867 

1.34 

-3.69 

-.03 

-2.00 

.781 

1 .95 

-1.17 

-.03 

-2.01 

.829 

1.50 

-3. 64 

-.04 

-2.49 

.71 1 

2.59 

-1.12 

-.05 

-2.51 

.794 

1.69 

-3.57 

-.07 

-2.99 

.642 

2.86 

-1.06 

-.08 

-3.0? 

.760 

1 .95 

-3.40 

-.11 

-3.49 

• 5R3 

3.15 

-.99 

-.10 

-3.52 

.724 

2.32 

-3.37 

-.15 

-3.98 

.532 

3.46 

-.90 

-.14 

-4.03 

.685 

2.74 

-3.2? 

-.21 

-4.47 

.488 

3.79 

-.80 

-.18 

-4 . 54 

.646 

2.88 

-3.05 

-.28 

-4.96 

.449 

4.13 

-.68 

-.21 

-5.04 

.615 

3.01 

-2.04 

-.37 

-5.46 

.416 

4.48 

-.55 

-.10 

-5.55 

.587 

3.15 

-2.60 

-.48 

-5.96 

.387 

4.84 

-.40 

-.37 

-6.06 

.561 

3.29 

-2.32 

-.61 

-6.46 

.36? 

5.22 

-.25 

-.45 

-6.57 

.537 

3.44 

-2.02 

-.77 

-6.97 

.339 

5.62 

-.08 

-.55 

-7.08 

.515 

3.59 

-1.67 

-.95 

-7.47 

.319 

6.03 

.11 

-.65 

-7.60 

.494 

3.74 

-1.27 

-1.17 

-7,97 

.3"1 

6 .45 

.31 

-.77 

-8.12 

.475 

3.90 

-.83 

-1.41 

-8.47 

.284 

6.89 

.5? 

-.90 

-8.65 

.458 

4.06 

-.34 

-1.70 

-8.97 

.270 

7.34 

.75 

-1 .^5 

-9.18 

.441 

4.2? 

.22 

-2.02 

-9.47 

.257 

7.81 






PROPERTIES OF COALFSCED JFT * 

= 

.49 Y = 

-1.54 

7= -9. 

33 

ll/U JO = 

XCOORD 

YCOORD 

ZCOOPO 

U 1 

DT A 






.49 

-1.54 

-0.33 

1 .ono 

1 .00 






• 80 

-1.77 

-9.S2 

.975 

1.10 






1 .32 

-2.02 

-10.31 

.951 

1 .23 






1 .80 

-2.29 

-in. 80 

.926 

1 .40 






2.35 

-2.57 

-11.30 

.901 

1 .62 






2.96 

-2.89 

-11.79 

.877 

1.70 






3.66 

-3.23 

-1?.?8 

.885 

1.76 






4.44 

-3.61 

-12.77 

.835 

1 .83 






5.31 

-4.02 

-13.27 

.817 

1 .90 






6.28 

-4.46 

-13.76 

.801 

1.97 






7.37 

-4.95 

-10.25 

.786 

2.04 






8.59 

-5.48 

-10.75 

.772 

2.11 






9.94 

-6.06 

-IS. 24 

.759 

2.18 






11.44 

-6.70 

-IS. 73 

.748 

2.25 






13.10 

-7.40 

-16.22 

.737 

2.3? 






14.95 

-8.16 

-16.72 

.727 

2.39 






17.00 

-9.00 

-17.21 

.717 

2.47 






19.28 

-9.91 

-17.70 

.708 

2.54 






21.79 

-10.92 

-18.19 

.700 

2.62 






24.57 

-12.02 

-18.69 

.692 

2.69 






27.65 

-13.23 

-19.18 

.685 

2.7? 






31 .06 

-14.56 

-19.67 

.678 

2.85 







(c) Jet Centerlines (Final Iteration) 
FIGURE A3. (Continued) 
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(d) Induced Flow Properties (Initial Approximation) 
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(e) Induced Flow Properties (Final Iteration) 


FIGURE A3. (Concluded) 
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APPENDIX B 


MANUAL FOR COMPUTER PROGRAM 
Description 

The program, which is a modified version of the Jet Flow Field computer pro- 
gram developed by Northrop Corporation under AFFDL contract F33615-69-C-1602 
(reference 12), evaluates the induced velocity field due to one or two jets exhausting 
into an arbitrarily directed crossflow. 

The equations of motion governing the development of each jet are integrated nu- 
merically for the position of the jet centerline, the nondimensionalized mean jet speed 
and the nondimensionalized major diameter of the ellipse which represents the jet cross 
section in the mathematical model. The set of first order differential equations is in- 
tegrated by means of a fourth order Adams predictor/corrector routine with a Runge- 
Kutta starting solution. 

The induced velocity components due to each jet at a given point are then calcu- 
lated by replacing each jet with a representative singularity distribution of sinks and 
doublets along the jet centerline. The contributions to the induced velocity components 
from the singularity distribution are summed over the length of each jet centerline. 

The velocity components due to each of the singularity distributions are additive at 
every control point where induced velocities are to be evaluated. 

For the two-jet configuration, the distance between the jet centerlines is tested 
and when intersection of the two jets is indicated, a coalesced jet is established from 
continuity and momentum considerations. The coalesced jet is treated as another 
independent jet in the computations for the induced velocity field. 

Jet exit velocity stratification effects may be treated by utilizing the velocity 
ratio and the jet exit diameter for the equivalent ideal nozzle or, alternatively, by 
accounting for the internal mixing which takes place through the introduction of input 
parameters A and B which are described in the discussion of the input data. These 
approaches to treating different types of jet nozzles producing stratified exit flows 
are discussed in detail in the section of the report dealing with velocity stratification 
effects. 
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For two-jet configurations, an iterative procedure involving modification of the 
crossflow into which the jets exhaust has been incorporated to account for further 
mutual interference effects between the two jets. The centerlines of the jets and their 
representative singularity distributions are calculated, using the unperturbed uniform 
crossflow in the computations. The jets are then broken into segments and the induced 
velocity, due to the other jet, is computed for each segment. The crossflow velocity 
vector for each segment is then modified by this induced velocity and the computations 
are repeated. In this iteration, and subsequently, each segment is treated as a sepa- 
rate jet exhausting into a uniform crossflow. 

Restrictions 

Jets must exhaust at some angle into the crossflow, i. e. the jet exhaust direction 
may not coincide with the freestream direction. 

For a two-jet configuration, the jet exits must both lie in the same XY plane, 
and the jet exhaust planes, defined by the freestream vector and the initial jet exhaust 
vectors, must be parallel (see figure B1 for definition of coordinate system). 

Control points at which the jet-induced velocity components are to be evaluated 
may not lie within the jet exhaust itself, as the formulation of the mathematical model 
is not valid in this region. 


Options 

• Induced Velocity Computation: Coordinates of the points at which velocity compo- 
nents are to be evaluated are provided as part of the input to the program. Only 
the induced velocities are computed at each point specified. 

• Induced Pressure Computation: Coordinates of the points at which the induced 
pressures are to be evaluated are provided as part of the input to the program. 

In adddition to the induced velocity components, the induced pressure in form 
of the flat plate pressure coefficient is evaluated at each point specified. 

• Note: If it is desired to use this modified version of the Jet Flow Field program 
in conjunction with the Transformation Method program described in Vol III of 
reference 12, some minor changes will be required. These consist primarily 
of including subroutines TRWING, TRBODY and ADAPT as part of the program 
and providing the input cards of Group B or Group C as described in Section II, 
Vol. Ill, reference 12. The program in its present form may be used to generate 
input to the Lifting Surface program described in Vol. Ill, reference 12 by 
exercising the punch control option. 
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Operating Information 


Core and Time Requirements: 

Computer: CDC 6600 

Core: 100 Kg to load 

61. 2 Kg to execute 

Time: Approximately one minute for a typical run using 250 control 

points and two iterations 

Additional Requirements: None 


Input Data 

The input data cards required by the program are shown in figure B2. The cards 
of Group I describe the jet configuration and provide parameters needed for computa- 
tional purposes. The cards of Group n describe the control points at which the jet- 
induced flow field is to be evaluated. 


Group I: Description of jet configuration, computational parameters 


Card 

No. 


© 


Variable 

Format 

Description 


MULT 

16 

Specifies number of jets in configuration 
MULT = 1, 2 

IGE0M 

16 

Specifies option of program being exercised 



1 

If IGE0M < 

1 

= 3 only induced velocities are 
evaluated 

= 4 flat plate pressure coeffi- 
cients are also evaluated 

IPUNCH 

16 

Punch control 



If IPUNCH \ 

( = 0 no punched output 

(. = 1 punched output generated 

NSEG 

16 

Number of integration intervals per jet segment 
Limit: 3^ NSEG?? 10 

ITER 

16 

Number of iterations to be performed on 


mutual interference velocities. May be left 
blank for single-jet configurations or two-jet 
configurations with the jet exits aligned in the 
freestream direction. For other two- jet con- 
figurations, ITER = 2 will normally suffice 
(see discussion, page 25) 
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Card 

No. 


© 


© 


© 


© 


© 


♦ 




Variable 

Format 

Description 

ALFA 

F12.0 

Angle of attack a (see figure Bl) | 

Angle of sideslip (3 (see figure Bl) j in degrees 

BETA 

F12.0 

N 

16 

Total number of intervals to be used in the nu- 
merical integration of the jet centerline 
Limit: N^lOO 

G 

F12.0 

Interval size to be used in the numerical inte- 
gration of the jet centerline, given as a fraction 
of the leading jet exit diameter. The integration 
routine will optimize the actual step size being 

— 


utilized. 

XJET 

F12.0 

X-coordinate of center of jet exit 

YJET 

F12.0 

Y-coordinate of center of jet exit 

ZJET 

F12.0 

Z-coordinate of center of jet exit 

PHI 

F12.0 

Jet exhaust angle 0 (see figure Bl) ) . , 

in degrees 

Jet exhaust angle <A (see figure Bl) ) 

PSI 

F12.0 

DJET 

F12.0 

Jet exit diameter 

^VELJ 

F12.0 

Freestream to jet exhaust velocity ratio 

A 

F12.0 

Ratio of effective core diameter to jet exit 
diameter for annular or vaned nozzles 
(see discussion, page 11, for details) 

B 

F12.0 

Jet mixing parameter for annular or vaned 

— m 


nozzles (see page 11) 


A and B must be set to zero for a nozzle with uniform exit flow or when treating 
stratification effects by using an equivalent ideal nozzle. A corresponds to the 
parameter a, and B corresponds to the parameter b in the discussion of annular 
and vaned nozzles (pp 10-14), 

Cards 4, 5, 6 are repeated to describe the second jet, if MULT = 2. For two- 
jet configurations, the upstream jet is listed ahead of the downstream jet. 
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Card 

No. Variable Format Description 


DIA 


<Z> 


F12.0 Empirical factor for coalesced jet. Function of 

jet orientation angle Q, which is the included 
angle between line connecting jet exit centers and 
the freestream velocity vector. If S2 < 20°, set 
DIA=1. 0. If S2220°, set DIA=0. 5. 

May be left blank for single jet configuration. 


Group II: Description of points where induced velocities/and pressures are computed 


(D 


[ 




NS 

16 

Number of spanwise control stations 

NC 

16 

Number of control points at each station 

X0(I) 

F12.0 

X-ccordinate of control point 1 

1 1 = 1, NCxNS 

Y0(I) 

F12.0 

Y-coordinate of control point ) Limit- 1^600 

Z0(I) 

F12.0 

Z-coordinate of control point ' 


Output 

Both printed and punched output may be obtained. 


Printed Output 

The jet configuration being treated is identified both by appropriate heading 
and by printout of pertinent input information. Jet centerline data printed out for 
all the jets in the configuration, including the jet resulting from the coalescence of 
two exhausting jets, consists of the centerline coordinates, the nondimensionalized 
mean jet speed, and the nondimensionalized major diameter of the ellipse representing 
the jet cross section. The point of intersection of the two exhausting jets in a two-jet 
configuration is identified, and the initial conditions for the resulting merged jet 
are given. 

The induced velocity components U, V, W, all nondimensionalized by U M , are 
printed out for each control point specified as part of the input. Additionally, if 
IGEOM = 4 was specified, the flat plate pressure coefficient, computed by using an 
image system, is printed out at each control point. 

For a two-jet configuration with the jets not aligned in the freestream direction, 
where a number of iterations are specified to account for mutual interference effects 
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between the two jets, the information described above is printed out for each iteration. 
An example of the printout for a typical problem involving the iterative process to 
account for mutual interference effects may be found in Appendix A. 


Punched Output 

Punched cards may be generated which can be utilized as part of the input to the 
Lifting Surface program described in reference 12, Vol in. The nondimensionalized 
velocity component W is punched out for every control point. This can serve as an 
approximation to the tangent of the jet-induced downwash angle for small angles of 
attack. Thus, the punched output from this option can serve as the downwash matrix 
[W] in the input to the Lifting Surface program. 


Programming Information 

Logical Structure 

The logical flow chart for the program is shown in figure B3. 
Purpose of Subroutines 


BITEST 

SEGMNT 

INTEG 

M0DIFY 

C0MP 

MUINT 

BALANC 

FIX 

0UTPT \ 
0UTPT1 f 

VEL0C 

DERIV 

PRT0UT 

TRANS1 1 
TRANS2 ( 

VEL1 


Tests for blockage and intersection of jets for two-jet configurations 
Breaks jet into segments for inclusion of mutual interference effects 
Integrates equations of motion for the jet path 
Computes mutually induced velocities 

Computes extent of overlap between jets in a two-jet configuration 
Computes modified freestream vector 

Establishes initial conditions for the coalesced jet from a momentum 
balance 

Limits maximum value of mutually induced velocities 

Transforms local coordinates to program coordinates 

Evaluates induced velocities at one control point due to one 
singularity distribution 

Computes derivatives for ADAMS 

Prints out computed answers 

Transforms input coordinates to program coordinates 

Computes effective velocity ratio for the downstream jet in a two- 
jet configuration 
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TRANS3 

PLANE 

ADAMS 

CFCAL \ 
CFCAL1 J 

R0TATE 

XPR0D 

S0L 


Transforms program coordinates to output coordinates 

Computes point of intersection between a given plane and a 
given line 

Adams predictor/corrector routine 

Computes direction cosines for jet-centered coordinate system 

Transforms program coordinates to jet-centered coordinates 

Computes cross product of two vectors 

Solves a system of three simultaneous equations 


Interdependence of Subroutines 


The Calling-Called matrix for the program is shown in figure B4. 





70 


FIGURE B2. INPUT DATA 




KT-KT +1 


MAIN PROGRAM 



FIGURE B3. LOGICAL FLOW CHART 
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M0DIFY 


VEL0C 


0UTPT 


0UTPT1 


ADAMS 


CFCAL 


FIGURE B4. CALLING-CALLED MATRIX 
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APPENDIX C 


LISTING OF COMPUTER PROGRAM 


PROGRAM STRJET 1 I NPUT , OUT PUT , PUNCH , T APE5» INPUT, TAPE6=0U TPUT , 

1 TAPE7=PUNCH) 

EVALUATION OF JET-INDUCED VELOCITY FIELD IMAXIMUM OF 2 JETS) 

BOTH JETS CAN HAVE STRATIFIED EXIT FLOW 

INITIAL JET EXHAUST DIRECTION MUST BE THE SAME FOR BOTH JETS 

DIMENSION XI ( 11,10) ,Z1UI,10) ,UJ1( 11 , 10 ) , D1 ( 11, 10 ) , DXDZ1I 11, 10 ) 
DIMENSION X2( l l, 10) ,22(11,10) ,UJ2( 11,10 ),D2( 11,10 ),DXDZ2( 11,10) 
DIMENSION X1T(100),Z1T( 100 )«UJ1T(100)«D1T(100),DXDZ1T( 100) 
DIMENSION X2T(100) ,Z2T( 100), UJ2TU00), D2T 1100), DXDZ2T1 100) 
DIMENSION X3( 100) ,Z3 (100) ,UJ3( 100) ,03(100 ),0XDZ3( 100) 

DIMENSION XBAS1 (100) , YBAS1 ( 100 ) ,ZBAS1 ( 100 ) 

DIMENSION XBAS2 ( 100 ) , YBAS2 ( 100) , ZBAS2 ( 100 ) 

DIMENSION XB SIT (100) , YBS1T( 100) , ZBSITl 100 ) 

DIMENSION XBS2T(100),YBS2T(100),ZBS2T(100) 

DIMENSION XBAS3 (100) , YBAS3 ( 100 ) , ZBAS3 ( 100 ) 

DIMENSION CF1(3,3,10),CF2(3,3,10),CF1T(3, 3),CF2T(3, 3),CF3( 3,3) 
DIMENSION UUEll 11,10) ,UUE2( 11,10 ),UUE1T (100) ,UUE2T( 100), UUE3( 100) 
DIMENSION SDXDZ11 11,10) ,SDXDZ2( 11, 10),SXZ1T( 100), SXZ2T( 100), 

1 S0XDZ3( 100) 

DIMENSION PARI 15) 

DIMENSION XJK10) ,YJ1(10) ,ZJ1 ( 10 ) , DJET1 ( 10 ) , VEL J 1 ( 10 ) 

DIMENSION XJ2I10) ,YJ2(10) ,ZJ2 ( 10 ) , DJET2 ( 10 ) ,VELJ2( 10) 

DIMENSION ALFQ1 ( 10 ) , BETQ1 ( 10 ) *GETQ1 ( 10 ) , ALFQ2 ( 10) , BETQ2 ( 10 ) , 

1 GETQ2 ( 10 ) 

DIMENSION DIRK 10) ,DIR2 ( 10) , ZS01 ( 10 ) ,ZS02 ( 10 ) , UFACTl ( 10 ) , 

1 UFACT2 (10) 

DIMENSION Ul( 10), VI (10) • W1 ( 10 ) ,U2 ( 10 ) , V2( 10 ) , W2 ( 10 ) 

C0MM0N/BLK1/CF1 ,CF2,CF1T,CF2T,CF3,UUE1,UUE2,UUE1T,UUE2T,UUE3,PAR 
C0MM0N/BLK2/X1 »Zl ,UJ1,D1,DXDZ1,X2, Z2,UJ2, D2,DXDZ2 
C0MM0N/BLK3/X1T ,Z1T,UJ1T,D1T,DXDZ1T*X2T,Z2T,UJ2T,D2T,DXDZ2T 
C0MM0N/BLKA/X3,Z3,UJ3,D3,DXDZ3 

C0MM0N/BLK5/XBAS1 ,YBAS1,ZBAS1 ,XBAS2,YBAS2, ZBAS2, XBAS3,YBAS3,ZBAS3 
COMMON /BLK6/XBS1T ,YBS1T,ZBSIT ,XBS2T ,YBS2T, ZBS2T 
COMMON/ 8LK7/ALFQ* BE TQ»GETQ.F1 ,F2,F3, VKONST 
COMMON /BLK8/ALF Ql ,BETQ1 ,GETQ1 , ALFQ2,BETQ2, GETQ2 
C0MM0N/BLK9/MULT , I HOLD 1 ,KOUNT 1 , I ONE, I TWO, ITHR,N1, N2, N3, IF I XI 
C0MM0N/BLK10/XJ1 ,Y J1 , ZJ1,DJET1,VELJ1,XJ2,YJ2,ZJ2,DJET2,VELJ2 
COMMON /BLKli /X JIT ,YJ1T,ZJ1T,DJET1T,VELJ1T,XJ2T,YJ2T,ZJ2T,DJET2T, 

1 VELJ2T 

COMMON /BLK12 /XJ3,YJ3 , Z J3,0JET3, VELJ3 
CQMM0N/BLK13/G «G2 ,G3,STEPI , STEP 1 2, STEP 1 3 
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C0MM0N/BLK14/V2X1 * V2YI » V2Zl»V2X2»V2Y2»V2Z2 

C0MM0N/BLK15/DI Rl *DIR2»DI R1T * DIR2T , DR3, ZS01* ZS02, UFACT1 ,UFACT2 
C0MM0N/8LK16/SDXDZ1 *SDXDZ2 *SDXDZ3* SXZ1T*SXZ2T 
C0MM0N/BLK17/GS*A,B,C,ISl,IS2,NPS 
C0MM0N/BLK18/U1 ,Vl*Wl,U2,V2*W2 
C0MM0N/BLK19/DI ARAT *DREF 
C0MM0N/BLK20/A1 ,B1 ♦ A2 ,B2 , ZSTORi , ZSTOR2 
C 

DIMENSION XO I 600) ,Y0(600) * ZO ( 600 ) » U( 600 ) * V ( 600 ) * Wl 600 ) 

DIMENSION CP ( 600) 

DIMENSION PH 10(3) » PS I 0(3) 

SET PARAMETERS 

El = .45 
E2 « .08 
E3 = 30. 

PI = 3.1416 
Cl = 2.24 

READ IN JET DATA 

READ (5*501) MULT * IGEOM* I PUNCH* NPS *NOI T 
READ (5*502) ALFA, BETA 
READ (5*503) N,GS 

501 FORMAT (1216) 

502 FORMAT (6F12.0) 

503 FORMAT (I6*F12.0) 

READ (5*502) X J1 ( 1 ) *Y J1 ( 1 ) ,ZJ 1 ( 1 ) , PH ID( 1 ) , PS ID( 1 ) , DJET 1 ( 1 ) , 

1 VELJl(l) 

READ (5,502) Al.Bl 
IF (MULT-2) 4,2,2 

2 READ (5,502) X J2 ( 1 ) , Y J2 ( 1 ) *ZJ2 ( 1 ) * PHID( 2 ) * PS I D ( 2 ) * DJET2 ( 1 ) * 

1 VELJ211) 

READ (5,502) A2*B2 
4 CONTINUE 

READ (5,502) DIARAT 
WRITE (6*690) 

690 FORMAT (1HI) 

IF (MULT-2) 14,15,15 

14 WRITE (6,603) 

603 FORMAT ( IH0,44X ,32H*** SINGLE JET CONFIGURATION ***/ > 

N 1 = N+l 

GO TO 17 

15 WRITE (6,604) 

604 FORMAT ( 1H0,45X ,29H*** TWO-JET CONFIGURATION ***/) 

17 CONTINUE 

WRITE (6,606) X Jl( l ) , YJ1 ( l ) , Z Jl ( l ) , PHID (1) , PS ID( 1 ) , VEL J l( 1 ) , 

1 DJETitl) 

606 FORMAT ( IHO, 15X.4HX JET, 11 X,4HYJET , 1 IX , 4HZ JET, 12X, 3HPHI , 12X, 3HPSI , 
112X,5HU/UJ0,llX,2HD0/8X,Fi5.4,lX,Fl4.4,lX,F14.4,lX,F14.4,lX,Fi4.4, 
21X,F14.4*1X,F14.4) 

IF (MULT-2) 20,18,18 

18 WRITE (6,607) X J2 ( 1 ) , YJ2 < 1 ) , Z J2 ( 1 ) , PHID(2 ) , PS ID ( 2 ) , VEL J2( 1 ) , 
l D JET2 ( 1 ) 
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607 FORMAT (8X,F15«4,1X,F14.4, 1X,F14.4, IX,F14.4, 1X,F14.4, IX, F 14. A, IX, 

1 F14.4) 

20 CONTINUE 

IF (Al) 41,41,42 

41 IF (MULT-2) 43,44,44 

43 WRITE (6,610) 

60 TO 56 

610 FORMAT ( 1H0, /15X,44HN02ZLE HAS UNIFORM EXIT FLOW CHARACTERISTICS) 

44 IF ( A2) 45,45,46 

45 WRITE (6,611) 

60 TO 56 

611 FORMAT ( 1H0, /15X, 51HB0TH NOZZLES HAVE UNIFORM EXIT FLOW CHARACTERI 
1STICS ) 

46 WRITE (6,612) A2,B2 
60 TO 56 

612 FORMAT ( 1H0, /15X,53HN0ZZLE OF JET 1 HAS UNIFORM EXIT FLOW CHARACTE 
IR I STICS/15X, 52HN0NUNIF0RM EXIT FLOW PARAMETERS FOR NOZZLE OF JET 2 
2: ,3X»4HA2 *,F5. 3,3X ,4HB2 =,F6.3) 

42 IF (MULT-2) 47,48,48 

47 WRITE (6,613) A1,B1 
60 TO 56 

613 FORMAT ( 1H0, /15X,47HN0NUNIF0RM EXIT FLOW PARAMETERS FOR THE NOZZLE 
1:,3X,4HA1 =,F5.3,3X,4HBl =,F6.3) 

48 IF (A2) 49,49,55 

49 WRITE (6,614) Al,81 
60 TO 56 

614 FORMAT ( 1H0,/15X,52HN0NUNIF0RM EXIT FLOW PARAMETERS FOR NOZZLE OF 

1 JET 1 : ,3X,4HAL =,F5.3,3X,4HB1 =»F6.3/15X, 53HN0ZZLE OF JET 2 HAS UN 
21 FORM EXIT FLOW CHARACTERISTICS) 

55 WRITE (6,615) A1,B1,A2,B2 

615 FORMAT ( IHO,/15X,52HNONUNIFORM EXIT FLOW PARAMETERS FOR NOZZLE OF 

1 JET 1 : ,3X,4HA1 =,F5.3,3X,4HB1 =,F6.3/15X, 52HN0NUNIFQRM EXIT FLOW P 
2ARAMETERS FOR NOZZLE OF JET 2:*3X,4HA2 =, F5 . 3, 3X, 4HB2 =,F6.3) 

56 CONTINUE 

WRITE (6,608) ALFA, BETA 

608 FORMAT ( IH0,/15X ,19HAN6LE OF ATTACK = , IX, F7.2/ 15X, 19HAN6LE OF SID 
IE SLIP =,1X,F7.2) 

WRITE (6,609) N,6S 

609 FORMAT { 1H0, / 15X ,32HNUMBER OF STEPS IN INTE6RAT ION =, IX, 13, /15X, 22H 
IINTE6RATI0N INTERVAL =,1X,F5.2,1X, 18HJET EXIT DIAMETERS) 

IF (MULT-2) 58,57,57 

57 WRITE (6,616) DIARAT 

616 FORMAT ( 1H0, / 15 X,88H INITIAL RATIO OF MINOR TO MAJOR DIAMETER OF EL 
1LIPTICAL CROSS SECTION OF COALESCED JET IS,F4.1) 

58 CONTINUE 

CALL TRANS1 ( MULT , ALFA, BETA, PS I D ) 

DO 31 1=1,10 
U 1 ( I ) = 0. 

VIII) = 0. 

Will) = 0. 

U 2 ( I ) = 0. 

V2(I) = 0. 

31 W2( I ) = 0. 

ITER = 0 
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*** START OF ITERATION LOOP **♦ 

30 CONTINUE 

DO 8 I«l,HULT 
PHI * PHI 01 1 )*• 0174533 
PSI = PSl0(I)*a0174533 
IF C 1-2) 5,6,6 

5 CONTINUE 

V2XI * SIN(PHI)*COS(PSI) 

V2Y1 = COSIPHI ) 

V2ZI = SIN(PHI)*SIN(PSI) 

CALL MUINT I ALFQ,BETQ,GETQ,U1 ( U , VI (1) , ML CD , ALFQ1 (1) , BETQ1 ( 1 ) , 
1 GETQ1( l) ,UFACT1( 1) ) 

CALL CFCALl ( ALFQ1 ( 1 ) ,BETQ1 I 1 ) ,GETQ1 t 1 ) , V2X1 , V2Y1, V2Z 1, CF1 , l ) 
CALL ROTATE ( V2XI , V2Y1, V2Z1,CF1 I 1, 1, 1 ) , VXT, VYT, VZT, 0 ) 

UJ1(1,1) - 1. 

01(1,1) * 1. 

Xl(l,l) * 0. 

Z 1 ( 1 • 1 ) « 0. 

0XDZ1 (1,1) * VXT/VZT 
XBASl(l) * XJ1 ( 1 ) 

YBASl(i) * YJ1 ( 1 ) 
zBAsim * zjim 
IF (ITER) 34,33,34 

33 A * CF1( 3,1, 1 ) 

B * CF1 ( 3 ,2, 1 ) 

C * CF1 ( 3 ,3, 1 ) 

34 CONTINUE 

COSTHP = 1* /SORT (!•+( VXT/VZT)** 2) 

COSTH *= A*V2X1 +B*V2Y1 *C*V2Zl 
G = GS 

G = G*COSTHP /COSTH 
STEPI * • 2*G 
DIRK 1 ) = 1. 

ZSOl(l) * 0. 

Bl » Bl/COSTHP 
D = ATANI VXT/VZT) 

IF (VXT) SOI , 902,902 

901 FI = *3*C0S( 0) 

GO TO 903 

902 FI = • 3/C0S( D ) 

903 CONTINUE 
GO TO 8 

6 CONTINUE 

V2X2 = SIN( PH I ) *COS ( PSI > 

V2Y2 » COS(PHI) 

V2Z2 = SIN(PHI)*SIN( PSI ) 

CALL MUINT ( ALFQ, BETQ,GETQ,U2 ( 1 ) , V2( 1) , W2 I 1 ) , ALFQ2I 1 ) , BETQ2( l ) , 
1 GETQ2 ( 1 I «UFACT2 ( II) 

CALL CFCALl ( ALFQ2I 1 ) ,BETQ2 ( 1 ) , GETQ2 ( 1 I , V2X2, V2Y2, V2Z2, CF2, 1 ) 
CALL ROTATE ( V2X2 , V2Y2, V2Z2,CF2 ( 1 , 1, 1 ) , VXT, VYT, VZT, 0 ) 

UJ2I 1,1) = 1. 

02(1,1) » 1. 
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X2< i » 1 ) » o. 

Z2( l, l) = o. 

DXDZ2(l,l) * VXT/VZT 
X8AS2 ( 1 ) * XJ2 ( 1 ) 

YBAS2C1) * YJ2UJ 
ZBAS2 ( 1 ) * Z J2 ( 1 ) 

COSTHP * l./SQRTd.-MVXT/VZT) **2 ) 

COSTH = A*V2X2 +B4V2Y2 +C*V2Z2 
62 = GS*DJETl (l)/DJET2(l) 

G2 * G2*C0STHP /COSTH 
STEPI2 = .2*G2 
D IR2( 1 ) = 1. 

ZS02( l) * 0. 

B2 = B2/C0STHP 
D = ATANI VXT/VZT) 

IF ( VXT) 904,905,905 

904 F 2 = «3*C0S( D) 

GO TO 906 

905 F2 = .3/COSID) 

906 CONTINUE 

8 CONTINUE 

IF (ITER) 9,7,9 
7 CALL VELl ( MULT , ALFA, VK1 ) 

IF (MULT-2) 12,11,11 

11 CONTINUE 

COMPUTE INITIAL OVERLAP 

CALL XPROO (V2X1,V2Y1,V2Z1,ALFQ, BET Q, GETQ,CFNX, CFNY, CFNZ ) 
CALL XPROO (V2X2,V2Y2,V2Z2,ALFQ,BETQ,GETQ,XT2,YT2,ZT2) 

CALL PLANE ( CFNX,CFNY ,CFNZ,XJ1 ( 1 ) , YJ1 ( l ) , ZJ1 <1 ) , XT2, YT2, ZT2, 
1 X J2( 1 ) , Y J2( 1 ) , Z J2 ( l ) ,XI,YI,ZI) 

01 ST = SQRT( (XI-XJ211) )**2 +1 YI-Y J2 (1 ) ) **2 ♦ ( Z I-Z J2 ( 1 ) ) **2 ) 

R = DJET1 ( 1 ) *• 5-DI ST 
FACT * (l.0+R/(DJET2(l)*.5))4.5 
IF (FACT-1.) 3,10,10 
3 IF (FACT) 10,10,13 
13 TEST1 * 0JET1(1)*.5+DIST 
TEST2 * 0 JET2 ( l ) *• 5 
IF (TESTI-TEST2) 23,10,10 
23 FACT = DJET1 ( 1 ) /DJET2 ( 1 ) 

10 OVLP = l.-FACT 
GO TO 9 

12 OVLP = 0. 

9 CONTINUE 
PARI 1 ) = El 
PAR( 2 ) * E2 
PAR13) = E3 
PAR( 7 ) = PI 
PAR( 8 ) = Cl 
PARI 9 ) = 1. 

PAR 111) = Ai 
PAR( 12) * Bi 
PARI 13) = A2 
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PAR( 14 J * 82 
ZSTORl = 0. 

ZST0R2 = 0. 

N2 = 0 
N3 * 0 
I HOLD 1 = 0 
K0UNT1 * 0 
TNEG = 8ETQ*V2Y1 
DREF = DJETl(l) 

IFIX1 = 0 
KSEG = 1 

00 50 1 = 1, N 

1 ONE * I 
I TWO = I 
VKONST = VK1 

IF ( MULT-2) 24,22,22 

24 IS1 = I-< KSEG-1 ) *NPS 
GO TO 27 

22 IF (IHOLDi-l) 25,25,21 

TESTS FOR BLOCKAGE AND INTERSECTION, PART OF INTEGRATION LOOP 

25 CALL BITEST U, TNEG, KSEG) 

IF ( IFIX1 ) 27,27,26 

26 CALL SEGMNT (I, KSEG, 1) 

NL = IFI Xl-( KSEG— 1 ) *NPS 

21 CALL SEGMNT (I,KSEG,2) 

27 CONTINUE 

INTEGRATION OF THE EQUATIONS OF MOTION FOR THE JET PATH 

CALL INTEG ( I , TNEG, KSEG) 

IF (IH0L01-2) 28,50,50 

28 IF U-N) 29,50,50 

29 IF (I-KSEG*NPS) 50,40,40 
40 CALL SEGMNT (I, KSEG, 3) 

50 CONTINUE 

IF ( IFIX1 ) 51,51,52 

51 NL = IS1+1 

52 IF l MULT-2) 60,53,53 

53 IF (0VLP-.01) 60,60,54 

54 CONTINUE 

CALL MODIFY < KSEG, TNEG, NL) 

CALL FIX (U1,V1,W1,U2,V2,W2,KSEG) 

IF (ITER) 60,60,70 
60 CONTINUE 

READING IN CONTROL POINTS WHERE INDUCED VELOCITIES WILL BE COMPUTE 

READ (5,501) NSMAX,NC 
NK * NSMAX*NC 

READ 15,502) (XO(I),YO(I) , Z0( I ) , 1*1, NK) 

70 CONTINUE 

CALL TRANS2 (Y0,Z0,NK) 
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EVALUATE INDUCEO VELOCITIES AT EACH CONTROL POINT 

IF (MULT-2) 90*91,91 
9i KTR1 * 0 
KTR2 = 0 

IF ( TNEG ) 96,96,97 

96 KTR2 = KOUNTI 
GO TO 90 

97 KTR1 * KOUNTI 
90 CONTINUE 

IF ( NL-2 ) 71,71,72 

71 KSEG1 = KSEG-1 
GO TO 75 

72 KSEGl = KSEG 
75 CONTINUE 

00 80 J S 1 ,NK 
U(J) = 0, 

V(J) * 0. 

W(J) = 0. 

DO 80 I=1«KSEG1 
PARI 6 ) = VELJ1I I ) 

PARI 5 ) = FI 
PARI 9 ) = DIR1II) 

IF I MULT-2) 81,82,82 

81 IF I I-KSEG) 113,114,114 

82 IF II-l) 111,112,111 

112 NF * NPS+1-KTR1 
GO TO 83 

111 IF (I-KSEG) 113,114,114 

113 NF = NPS+i 
GO TO 83 

114 NF = NL 

83 CONTINUE 

CALL VELOC 1 1 ,NF,Z1 1 1 , 1 ) ,X1( 1 , I ) , DXDZ1 1 1, I ) ,UJ 1 1 1, I ) , 01 1 1, I ) , 

1 UUEl(l,n,XJl(I),YJl(I),ZJl( I),DJETlin,CFl(l,l,I),PAR,XOU), 

2 YOI J ) , ZOI J ) ,UI NO, VINO, WIN0,S0X0Z1 < 1, I ) ) 

UIJ) = U( J) +UIND 

VIJ) = VI J) +VIND 
WIJ) = WIJJ+WINO 
IF (MULT-2) 80,86,86 

86 CONTINUE 

PARI 6 ) = VELJ2II) 

PARI5) = F2 
PAR 1 9 ) = 0IR2II) 

IF (1-1) 87,115,87 

115 NF * NPS+1-KTR2 

87 CONTINUE 

CALL VELOC 1 1 ,NF, Z2 1 1 , 1 ) ,X2 (1 , I ) , DXDZ2U, I) ,UJ2 1 1, I ) , 021 1 , I ) , 

1 UUE2(1,I),XJ2(I),YJ2(I),ZJ2( I ) , DJET2I I), CF2 1 1, 1, I ) , PAR, XO I J ) , 

2 YOI J), ZOI J) ,UINO, VINO, WIND, SDXDZ2I 1,1) ) 

UIJ) * U! J) +UI NO 

VIJ) = V ( J) +VI NO 
WIJ) * W(J)+WIND 
80 CONTINUE 
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IF ( I H0LD1-1 ) 88*88,89 
89 CONTINUE 

PARI 6 ) = VELJ3 
PAR( 5 ) = F3 
PARI9 ) * 0R3 
N3 = ITHR+1 
00 120 J*1,NK 

CALL VELOC l 1 ,N3 , Z3 *X3,0XDZ3,UJ3, D3,UUE3, X J3, YJ3, Z J3, DJET3, 
1 CF3,PAR,X0(J) , Y0 I J ) , Z0 ( J ) *UIN0, VINO, WIND, SDX0Z3 ) 

UIJ) = U( JI+UIND 
V(J) = V( JI+VIND 
120 HIJ) * W( J)+WIND 
88 CONTINUE 

IF I IGEOM-3) 93,93,92 

COMPUTE FLAT PLATE PRESSURE COEFFICIENTS 

92 DO 85 J=1,NK 

CPT = 4 . *IUI J ) *1 ALFQ+UI J) ) +HI J )*| GETQ+W ( J )) ) 

85 CPU) = i.-(ALFQ*ALFQ +6ETQ*6ETQ +CPT ) 

93 CONTINUE 

CALL TRANS3 ( YO , ZO, V, W,NK, KSEG, NPS ,TNEG,NL ) 

PRINT OUT COMPUTED RESULTS 

CALL PRTOUT I IGEOM, XO, YO, ZO ,U, V,H,CP,NK, ITER,NOIT,OVLP ) 

IF I0VLP-.01) 94,94,98 
98 IF IOVLP-,99) 77,77,79 

77 DO 78 1=1 , KSEG 

uiin = ui(I)*ovlp 
vim = vi ( i ) *ovlp 
wim = wnn*ovLP 
U2II) = U21 I ) *0VLP 
V2 ( I ) = V2( I ) *0VLP 

78 W2(I) = W2 ( I ) *0VLP 

79 CONTINUE 

IF IITER-NOIT) 35,94,94 
35 ITER = ITER+l 
GO TO 30 

*** ENO OF ITERATION LOOP *** 

94 CONTINUE 

IF IIPUNCH) 95,99,95 

PUNCH OUT OATA FOR LIFTING SURFACE PROGRAM 

95 00 101 1*1, NK 
101 HU) * -HII) 

J1 = 1 
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00 102 I *1 *NSMAX 
J2 = JKNC-1 

WRITE (7,710) (W( J) ,J=J1, J2) 
102 J1 * J2+1 
710 FORMAT (5E14.7) 

99 CONTINUE 
STOP 
END 


SUBROUTINE SEGMNT <I,KSEG,IND) 


I ND=1 , ESTABLISHES INITIAL CONDITIONS 
IN0=2, INTEGRATES CONTINUATION JETS 
I ND=3 , ESTABLISHES INITIAL CONDITIONS 


FOR CONTINUATION JETS 
FOR A NEW SEGMENT 


C 


C 


EXTERNAL OERIV 


X1(11*10),ZK11«10),UJ 1(11,10), D1(L1,10)«DXDZK11«10) 
X2< 11,10) ,Z2( 11,10) ,UJ2 ( 11, 10 )« 02 ( 1 1, 10 ) , DX0Z2( 1 1, 10) 
X1T ( 100 ) , Z1T( 100),UJ1T(100),D1T(100), DXOZ 1T( 100) 
X2T(100),Z2T(100),UJ2T(100) , D2T (100),0XDZ2T(100) 

X3( 100) ,Z3(100) ,UJ3(100) ,03(100 ),DXDZ3( 100) 

XBASl(lOO) ,YBAS1( 100) ,ZBAS 1(100) 

XBAS2 ( 100 ) , YBAS2 (100), ZB AS 2(100) 

XBSIT(IOO) , YBS1T(100),ZBS1T(100) 

XB S2T ( 1 00 ) , YBS2T ( 100 ) ,ZBS2T(100 ) 

XBAS3 ( 100 ) , YBAS3 ( 100 ) , ZBAS3( 100 ) 

CF1 (3,3,10) ,CF2(3,3,10),CF1T(3, 3),CF2T(3, 3),CF3(3,3) 
UUE1( 11,10) ,UUE2( 11,10 ),UUE1T( 100) ,UUE2T( 100) ,UUE3( 100) 
SDXDZK 11,10) , SDXDZ2 ( 1 1 , 10 ) , SXZ IT ( 100 ) , SXZ 2T( 100 > , 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
1 SDXDZ3 (100) 

DIMENSION PARI 15) 

DIMENSION XJl(lO) , Y J1 ( 10) , ZJl ( 10 ) , DJET1 ( 10 ) , VEL J1 ( 10 ) 
DIMENSION XJ2(10),YJ2(10) , ZJ2 (10),DJET2(10),VELJ2(10) 
DIMENSION ALFQ1 (10) ,BETQ1 ( 10) ,GETQ1 ( 10) , ALFQ2( 10) • BETQ2 ( 10 ) , 
1 GETQ2 ( 10) 

DIMENSION DIRK 10) , DIR2 ( 10 ) , ZSOl ( 10) , ZS02 1 10 ) ,UFACT1 ( 10), 

1 UFACT2 ( 10) 

DIMENSION Ul( 10), VI (10) ,W1 ( 10 ) ,U2 1 10 ) , V2( 10) , W2 ( 10) 


C0MM0N/BLK1/CF1 ,CF2,CF1T,CF2T , CF3,UUE1,UUE2,UUE1T ,UUE2T ,UUE3,PAR 
COMMON/BLK2/X1 , Z1 ,UJ1 , D1 « DXDZ 1 ,X2« Z2,UJ2, D2, DXDZ2 
C0MM0N/BLK3/X1T ,Z1T,UJ1T,D1T,DXDZ1T, X2T , Z2T, U J2T, D2T , DXDZ2T 
C0MM0N/8LK4/X3, Z3 ,U J3 ,D3« DXDZ3 

C0MM0N/BLK5/XBAS1 , YBAS1 ♦ ZBAS1 , XBAS2, YBAS2,Z8AS2» XBAS3, YBAS3, ZBAS3 
C0MM0N/BLK6/XBS1T ,YBS1T,ZBS1T , XBS2T , YBS2T, ZBS2T 
C0MM0N/BLK7/ALFQ,BETQ«GETQ,F1 ,F2,F3, VKONST 
C0MM0N/BLK8/ALFQ1 ,BETQ1 ,GETQ1 , ALFQ2,BETQ2,GETQ2 
C0MM0N/BLK9/MULT , IH0LD1 ,K0UNT1, I ONE, ITWO, ITHR,N1,N2,N3, IFIX1 
COMMON/BLKIO/XJI , Y J1 , ZJl, 0JET1, VELJ1, XJ2, YJ2,ZJ2,DJET2,VELJ2 
C0MM0N/BLK11 /XJIT,YJ1T,ZJ1T,DJET1T,VELJ1T,XJ2T,YJ2T,ZJ2T,DJET2T, 

1 VELJ2T 

COMMON/BLK12/XJ3, YJ3 ,ZJ3,DJET3» VEL J3 
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C0MM0N/BLK13/G,G2,G3,STEPI .STEPI2, STEPI3 
C0MM0N/BLK14/ V2XI »V2Y1»V2Z1*V2X2»V2Y2,V2Z2 

COMMON/BLK15/OIRI ,DIR2» DIRIT * DIR2T ,DR3, ZS01, ZS02, UFACT1 » UFACT2 
C0MM0N/BLK16/SDXDZI *SDXDZ2,SDXDZ3, SXZIT.SXZ2T 
C0MM0N/8LK17/GS , A,B,C , IS1.IS2.NPS 
C0HM0N/BLK18/U1 ,V1 , Wi t U2 tV2 t W2 
COMMON/BLK20/A1,B1,A2,B2,ZSTOR1,ZSTOR2 
c 

DIMENSION FIN(4) ,F0UT(4) 

C 

IF CIND-2) 26,21.40 
26 NL = IFIXl-( KSEG— 1 ) *NPS 
X JIT = XBAS1 ( IONE ) 

YJ1T = YBAS1 < I ONE ) 

Z JIT = ZBAS1 ( I ONE ) 

UJ1T(1) = 1. 

DlT(i) = 1. 

X1T(1) = 0. 

ZlT(l) = 0. 

CALL CFCAL ( ALFQ,BETQ,GETQ,V2X1,V2Y1,V2Z1,CF1T ) 

CALL ROTATE ( V2X1 ,V2Y1 , V2Z1 ,CF1T, VXT, VYT, VZT, 0 ) 

DXDZlT(l) = VXT /VZT 
XBS1TC1) = XJIT 
YBSlT(l) = YJ1T 
ZBSlT(l) * Z JIT 

DJET1T = Dl( NL.KSEG) *D JET 1 (KSEG) 

VELJ1T = UJ1(NL,KSEG)*VELJ1(KSEG> 

GIT = GS*DJET1(1)/0JETIT 
STEPIT = GIT*. 2 

IF (DIRl (KSEG)— .2501) 320,320,321 

321 ZOVM * { ZS01 ( KSEG) +Z1 ( NL.KSEG) ) / ( VELJ1 ( KSEG)*UUE1 (NL-l, KSEG) ) 
IF ( Z0VM-F1 ) 314,314,320 
314 DIRIT = l 75*Z0VM/Fl 
GO TO 325 
320 DIRIT * .25 
325 CONTINUE 

XJ2T = XBAS2 ( I TWO) 

YJ2T * YBAS2 ( I TWO) 

ZJ2T = ZBAS2 ( ITWO) 

UJ2TJ1) = 1. 

D2T( 1 ) = l. 

X2T ( 1 ) * 0. 

Z2T( 1 ) = 0. 

CALL CFCAL ( ALFQ, BETQ.GETQ, V2X2, V2Y2, V2Z2.CF2T ) 

CALL ROTATE (V2X2,V2Y2,V2Z2,CF2T, VXT, VYT, VZT, 0) 

DXDZ2T ( 1 ) » VXT /VZT 
XBS2T ( l ) * XJ2T 
YBS2T ( l ) = YJ2T 
ZBS2T ( 1 ) a Z J2T 

DJET2T a 02 ( NL.KSEG) *DJET2 (KSEG) 

VELJ2T = UJ2 ( NL ,KSEG ) *VEL J2(KSEG) 

G2T = GS*DJET1 (D/DJET2T 
STEP2T a G2T*.2 
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IF I DIR2 I KSEG )-• 2501 ) 330,330,331 

331 ZOVM = ( ZS02 ( KSEG) +Z2I NL, KSEG) ) / I VELJ2 C KSEG ) *UUE2 INL-1,KSEG) ) 

IF { Z0VM-F2 ) 324,324,330 
324 DIR2T = 1 75*Z0VM/F2 
GO TO 335 
330 DIR2T * .25 
335 CONTINUE 
GO TO 50 
21 CONTINUE 

IA = I-IFIXl+1 
PAR ( 6 ) = VEL JIT 
PARI 5 ) = FI 
PAR ( 9 ) * DIR1T 
PARI 15) * 3. 

UUE1T ( IAI » 1. 

ZITlIA+ll * Z1T I IAJ+G1T 
F INI 1 ) = UJ1TIIA) 

F INI 2 ) * OiTIIA) 

F INI 3 ) * XiTIIA) 

FINI4) = 0XDZ1T 1 1 A) 

CALL ADAMS I4,Z1TCIA) , Z1T I IA+l ) , STEP IT, GIT, 999, 1.0E-04, 1.0E-05,0, 
1 FIN,FOUT,PAR,DERIV) 

UJ1TIIA+U = FOUTU) 

DlTIIA+1) * F0UTI2J 
XlTIIA+ll « F0UTI3) 

DXDZ1TIIA+1I = F0UTI4) 

SXZ1T I IA+l ) = PARI 10) 

CALL OUTPT IX1T(IA+1),Z1T(IA+1),DXDZ1T(IA+1),CFIT,0JET1T«XJ1T, 

1 YJIT,ZJIT,XBS1TIIA+1),YBS1TIIA+1),ZBS1T(IA+1), OUM, DUM, OUM ) 

PARI6) * VELJ2T 
PARI 5 ) = F2 
PAR 1 9 ) * DIR2T 
PARI 15) » 3. 

UUE2TIIA) = l. 

Z2TIIA+1) * Z2TI IAJ+G2T 
F INI 1 1 = UJ2TIIA) 

F INI 2 l = D2TIIA) 

F INI 3 ) = X2TI IA) 

F INI 4 ) = DX0Z2T I IA ) 

CALL ADAMS I4,Z2TCIA) ,Z2T I I A+l > ,STEP2T,G2T, 999, 1.0E-04, 1 .0E-05, 0, 
1 FIN,FOUT , PAR, DERI V) 

UJ2TI I A+l ) * FOUTU) 

D2TIIA+1) * F0UTC2) 

X2T I IA+1 ) « F0UTI3) 

DXDZ2T I IA+1 ) = F0UTI4) 

SXZ2TI IA + 1) = PARI 10) 

CALL OUTPT tX2T(IA+l),Z2T!IA+l), DXDZ2T IIA+i),CF2T,DJET2T,XJ2T, 

1 YJ2T*ZJ2T,XBS2T(IA+ll* YBS2T I IA+1),ZBS2TI I A+l ), DUM, OUM, DUM ) 

GO TO 50 
40 CONTINUE 

KSEG = KSEG+1 
X Jl( KSEG ) * XBAS1 1 1 ONE+1 ) 

YJ1IKSEG) * YBAS1 I IONE+1 ) 

Z J1IKSEG) * ZBASi I IONE+1 ) 


83 



UJKltKSEG) = 1. 

Old, KSEG) * 1. 

Xl( l, KSEG) = O. 

Z 1( l» KSEG) = 0. 

CALL MUINT ( ALFQ,BETQ,GETQ,U1 ( KSEG ) ,V1 (KSEG) ,W1 ( KSEG ), ALFQ 1 ( KSEG) , 
1 BETQl(KSEG) ,GETQ1( KSEG), UFACT1 t KSEG) ) 

CALL CFCAL1 ( ALFQ1 (KSEG) »BETQ1 (KSEG) »GETQ1 (KSEG) ,V2X1,V2Y1,V2Z1, 
l CFlt KSEG) 

CALL ROTATE ( V2XI,V2Y1,V2Z1,CF1 < 1, 1,KSEG) ,VXT,VYT, VZT, 0) 

DXDZ1 ( 1 *KSEG ) = VXT/VZT 

OJETUKSEG) = D1(IS1+1»KSEG-1)*DJET1 1 KSEG- 1 ) 

VELJ1 (KSEG) = UJKIS1 + 1 *KSEG-i ) *VELJ1 ( KSEG-1 ) 

COSTHP = 1./SQRTI1.-M VXT/VZT) **2 ) 

COSTH = A*V2X1 +B*V2Yl +C*V2Z1 
G * GS^OJETl ( 1 ) /0JET1 (KSEG ) 

G » G*COSTHP /COSTH 
STEPI = G*.2 

!F ( PAR ( 11 ) ) 423,423,422 

422 ZST0R1 = ZSTORl+ZM I Sl+1 » KSEG-1 ) *DJETi ( KSEG— 1 ) /DJET 1 ( 1 ) 

PARC 11) * ( A1-B1*ZST0RI )*DJET1 ( 1 ) /DJET 1 (KSEG) 

423 CONTINUE 

IF (DIRl(KSEG-l)-.250i ) 420,420,421 
421 ZOVM = ( ZS01 ( KSEG-1 )+Zl(ISi+l ,KSEG-1 ) ) / ( VELJ1 (KSEG-1 } * 

I UUEK ISi, KSEG-1 ) ) 

IF ( ZOVM— FI ) 414,414,420 
414 DIRK KSEG) = l.-,75*Z0VM/Fl 

Z SOI ( KSEG ) * (1«— DIR1 (KSEG) )*VELJ1 ( KSEG )*U FACT 1 (KSEG)* 

I UUEi (I Si, KSEG-1 )/UFACTl (KSEG-1 )*F1/. 75 

GO TO 425 

420 DIRK KSEG I = ,25 
425 CONTINUE 

IF (MULT-2) 50,450,450 
450 CONTINUE 

XJ2 ( KSEG ) = XBAS2 ( ITWO+1 I 
YJ2(KSEG) = YBAS2( ITWO+1) 

Z J2(KSEG) * ZBAS2 ( ITWO+1 ) 

UJ2I 1 » KSEG) = 1. 

02(1, KSEG) * 1. 

X2( 1»KSEG) * 0. 

Z2(i,KSEG) = 0. 

CALL MUINT ( ALFQ, BETQ,GETQ,U2 (KSEG ) ,V2( KS EG) ,W2( KSEG ), ALFQ2( KSEG ) , 
1 BETQ21KSEG) , GETQ2C KSEG ) ,UFACT2 ( KSEG) ) 

CALL CFCAL1 ( ALFQ2C KSEG) , BETQ2 (KSEG) ,GETQ2 (KSEG) , V2X2, V2Y2,V2Z2, 

1 CF2.KSEG) 

CALL ROTATE ( V2X2 ,V2Y2, V2Z2.CF2 ( 1, 1, KSEG) , VXT,VYT, VZT, 0 ) 

DXDZ2 ( 1 , KSEG ) * VXT/VZT 

D JET2 ( KSEG) = D2( I S2+1, KSEG-1 ) *DJET2 (KSEG-1 ) 

VEL J2 (KSEG) * UJ2( IS2+1 , KSEG-1) *VELJ2(KSEG-i ) 

COSTHP = l./SQRT(l»+( VXT/VZT) **2 ) 

COSTH = A*V2X2 +B4V2Y2 *C»V2Z2 
G2 = GS*0JET1 ( 1 ) /DJET2( KSEG) 

G2 = G2*COSTHP/COSTH 
STEPI 2 = G2*.2 
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non 


IF ( PARI 131 ) 433*433,432 

432 ZST0R2 = ZST0R2+Z2 II S2+1 * KSEG-i ) *DJET2( KSEG-l ) /0JET2 { I ) 
PARI 13) * I A2-B2*ZST0R2 ) *DJET2 ( 1)/DJET2(KSEG) 

433 CONTINUE 

IF (0IR2(KSEG-1)-.2501) 430,430,431 
431 ZOVM = ( ZS02 ( KSEG-l ) +Z2 1 1 S2+1 , KSEG-l ))/(VELJ2(KSEG-l)* 

1 UUE2I I S2» KSEG-l ) ) 

IF (Z0VM-F2) 424,424,430 
424 0IR2IKSEG) = l.-.75*Z0VM/F2 

ZS02IKSEG) * < 1 .-01 R2 ( KSEG ) ) *VELJ2 ( KSEG ) *UFACT2 ( KSEG ) * 

1 UUE2IIS2,KSEG-I) /UFACT2 (KSEG-i )*F2/. 75 

GO TO 435 

430 0IR2IKSEG) * .25 
435 CONTINUE 
50 CONTINUE 
RETURN 
END 


SUBROUTINE MODIFY ( KSEG.TNEG, NL ) 

COMPUTES MUTUALLLY INDUCED VELOCITIES 

DIMENSION X1(11,10)*Z1(11, 10 ),UJ1(11«10), 01(11,10), DXDZ 11 11,10) 
DIMENSION X21 11,10) ,Z2( 11,10) ,U J2 1 11, 10 ) , D2 1 11, 10 ) • DXDZ21 1 1, 10 ) 
DIMENSION XIT(IOO) ,Z1T( 100),UJ1T( 100),D1T 1 100) , DXDZ IT 1 100 ) 
DIMENSION X2TU00) , Z2T 1 100) ,U J2T 1 100 ) ,D2T ( 100) , DXDZ2T 1 100 ) 
DIMENSION X3I 100) ,Z3(100) ,UJ3< 100) ,03( 100 ) ,DXDZ3( 100) 

DIMENSION XBAS1 (100) ,YBAS1 1100) , ZB ASH 100) 

DIMENSION XBAS2 I 100 ) , YBAS2 ( 100) , ZBAS2I 100 ) 

DIMENSION XBSIT(IOO) ,YBS1T(100) ,ZBS1T(100) 

DIMENSION XBS2TI100) ,YBS2T(100) ,ZBS2T( 100 ) 

DIMENSION XBAS3I100) ,YBAS3I100) ,ZBAS3( 100 ) 

DIMENSION CF1(3,3,10),CF2(3,3*10),CF1T<3,3 ) ,CF2T (3,3),CF3I3,3) 
DIMENSION UUE 1(11,10) ,UUE2 (11,10),UUE1T (100) «UUE2T 1 100) ,UUE3( 100) 
DIMENSION SDXDZ1I 11,10) *S0XDZ2( 11, 10 ) , SXZ IT ( 100 ) , SXZ2TI 100), 

1 SDXDZ3( 100) 

DIMENSION PAR( 15) 

DIMENSION XJl(lO) ,YJ1(10) ,ZJ1 ( 10 ) , DJET1 ( 10 ) , VEL J 1 ( 10) 

DIMENSION XJ2 (10) , Y J2 ( 10) , ZJ2 I 10 ) , DJET2 ( 10 ) , VEL J2( 10) 

DIMENSION ALFQl(lO), BET 01(10) ,GETQ1 ( 10 ) , ALFQ2I 10 ) ,BETQ2 ( 10 ) , 
l GETQ2I 10) 

DIMENSION DIRK 10) ,DIR2 ( 10) ,ZS01 ( 10) ,ZS02 ( 10),UFACT1( 10), 

1 UFACT2 ( 10) 

DIMENSION Ul( 10) ,V1(10) ,H1(10),U2( 10),V2( i0),W2(10) 

COMMON /BLKl/CFl ,CF2,CF1T,CF2T,CF3,UUE1, UUE2, UUEIT , UUE2T, UUE3,P AR 
C0MM0N/BLK2/X 1 , Z1 ,UJl ,D1 , DXDZ 1 , X2 , Z2 ,UJ2, D2, DXDZ2 
C0MM0N/BLK3/X1T »Z1T,UJ1T,D1T,DXDZ1T,X2T,Z2T,UJ2T«D2T«0XDZ2T 
COMMON/BLK4/X3, Z3,UJ3,D3,DXDZ3 

C0MM0N/BLK5/XBAS1 ,YBAS1,Z8AS1 ,X8AS2, YBAS2, ZBAS2, XBAS3,YBAS3,ZBAS3 
C0MM0N/BLK6/XBSIT »YB$1T»ZBS1T ,XBS2T, YBS2T , ZBS2T 
C0MM0N/BLX7/ALFQ,BETQ*GETQ,F1 ,F2,F3, VKONST 
C0MM0N/BLK8/ALFQ1 ,BETQ1 ,GETQ1 , ALFQ2,BETQ2,GETQ2 
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non 


C0MM0N/BLK9/MULT » IH0LD1 *K0UNT1 , I ONE* ITWO, ITHR«Nl*N2*N3* IFIX1 
C0MM0N/BLK10/XJI, YJ1* ZJ1 ,D JET1 , VEL J1*XJ2, YJ2, ZJ2, DJET2, VEL J2 
COMMON/BLKll /XJ1T *YJ1T*ZJ1T*DJET1T » VEL JIT » XJ2T* YJ2T* ZJ2T* DJET2T* 
1 VELJ2T 

COMMON/BLK12/X J3, YJ3* ZJ3,DJET3* VEL J3 
C0MM0N/BLK13/G,G2,G3,STEPI,STEPI2,STEPI3 
C0MM0N/BLKIA/V2X1 * V2Y1 , V2Z1 * V2X2* V2Y2* V2Z2 
C0MM0N/BLKi5/DIRl,0IR2,DIRiT,DIR2T,DR3,ZS01*ZS02,UFACTl,UFACT2 
C0MM0N/BLK16/SDXDZ1 »SDXDZ2*SDXDZ3*SXZ1T*SXZ2T 
C0MM0N/8LK17/GS * A »B *C » ISl*IS2tNPS 
COMMON/BLK18/Ul,Vl,Wl,U2,V2,W2 
C 

DIMENSION XPltlO)*YPH10)*ZPH10)»XP2I10),YP2{10)*ZP2(10) 

CHOOSING POINTS ON CENTERLINE 

KTR1 = 0 
KTR2 = 0 

IF ITNEG1 206*206*207 

206 KTR2 = KOUNTI 
GO TO 208 

207 KTR1 = KOUNTI 

208 CONTINUE 

DO 70 I=l,KSEG 
IF (I-1J 71*71,72 

71 IK = l 
GO TO 73 

72 IK = I I-l)*NPS*l-KTRl 

73 CONTINUE 

XPlt I ) = XBAS1I IK) 

YPl ( I ) = YBASH IK) 

ZPl(I) = ZBASK IK) 

U1U) = 0. 

VUI) = 0. 

W1U) = 0. 

70 CONTINUE 

DO 80 I*1,KSEG 
IF (1-1) 81,81,82 

81 IK * 1 
GO TO 83 

82 IK = I l-i)*NPS+l-KTR2 

83 CONTINUE 

XP2III = XBAS21 IK) 

YP2II) = YBAS2I IK) 

ZP21I) = ZBAS21IK) 

U2(I) = 0. 

V 2 < l ) » 0* 

W2II) = 0. 

80 CONTINUE 

PARIS) = FI 
DO 110 J s l *KSEG 
DO 110 I S 1 * KSEG 
PARI 6) = VELJ1II) 

PARI 9) = 0IR1II) 
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IF CI-1) 111*112*111 

112 NF 3 NPS+l-KTRl 
60 TO 115 

III IF CI-KSEG) 113*114,114 

113 NF * NPS+l 
60 TO 115 

114 NF = Nl 

IF (NL-2) 110,110*115 

115 CONTINUE 

CALL VELOC C 1 ,NF*Z1 C 1 , 1 ) *X1 <1 * I ) ,DXDZ1 C 1* I ) *UJ1 C 1, I) *Dl C 1 * I ) • 

1 UUEl(l,I),XJUI)*YJl(I),ZJim*DJETim*CFl(l,l*I)*PAR,XP2( J), 

2 YP2(J)*ZP2(J) *UIND* VINO* WIND* SDXDZ1 C 1,1)) 

U2CJ) = U2CJ1+UIND 

V2C J) 3 V21 J) +VIND 
W2CJ) = W2( J ) +WIND 
110 CONTINUE 

IF ( IH0LD1— 1 ) 150*150*120 
12C PAR ( 6 ) 3 VELJ3 
PAR ( 5 ) 3 F3 
PAR( 9 ) 3 DR3 
N3 3 ITHR+1 
DO 121 J s l*KSEG 

CALL VELOC C 1 *N3 » Z3 » X3* DXDZ3, UJ3* D3* UUE3* XJ3*YJ3*ZJ3*DJET3» 

1 CF3, PAR « XP2 C J ) *YP2 ( J) « ZP2CJ) *Uf NO* V IND* W IND* SDXDZ3 ) 

U2 ( J ) = U2C J ) +UI ND 
V2(J) 3 V2CJJ+VIND 
W2 ( J ) * W2C J)+WIND 

121 CONTINUE 

PAR ( 6 ) * VELJ2T 
PAR( 5 ) 3 F2 
PAR ( 9 ) 3 0IR2T 
DO 122 J 3 1*KSEG 

CALL VELOC (1 ,N3,Z2T,X2T,DXDZ2T,UJ2T,D2T*UUE2T,XJ2T, YJ2T,ZJ2T, 

1 DJET2T *CF2T*PAR*XP2 ( J ) *YP2( J ) » ZP2 C J ) , U IND* V IND* WIND* SXZ2T ) 
U2CJ) 3 U2( J) -UINO 
V2C J ) = V2(J)-VIN0 
W2 ( J ) 3 W2 ( J )-WI ND 

122 CONTINUE 
150 CONTINUE 

PAR{ 5 ) = E2 
DO 210 J 3 l,KSE6 
DO 210 I 3 1*KSEG 
PARC 6 ) 3 VELJ2CI) 

PARC 9 ) » DIR2C I ) 

IF (1-1) 211,212*211 

212 NF 3 NPS+1-KTR2 
GO TO 215 

211 IF CI-KSEG) 213*214*214 

213 NF 3 NPS+l 
60 TO 215 

214 NF 3 NL 

IF ( NL-2) 210*210,215 

215 CONTINUE 
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noon 


CALL VELOC <1 ,NF , 12 ( 1 • I ) , X2(l , I > # DXDZ2<1, I ) ,U J2(l, I > ,02 11 , I) , 

1 UUE2 (1,I),XJ2(I),YJ2(I), ZJ2 ( I ) , DJET2 ( I ) « CF2 (l«l,I),PAR,XPl(J)« 

2 YPK J) ,ZP1( J) ,UIND,VIND,WIND,SDXDZ2(1, I) ) 

Ul< J) = Ul( J) +UIND 

VKJ) = VK J ) 4-VINO 
W1(J) = WKJJ+WIND 
210 CONTINUE 

IF ( I HOLD l—l J 250,250*220 

220 PARI 6 ) * VELJ3 
PARI 5 ) = F3 
PARI9) = DR3 

00 221 J=1,KSEG 

CALL VELOC 1 1 ,N3 , Z3 ,X3, 0X0Z3,UJ3* 03,UUE3, XJ3, YJ3, Z J3, 0JET3, 
l CF3,PAR,XP1( J) * YPl I J ) « ZP1 IJ) ,UIND,VIND,WIND,SDXDZ3) 

UltJ) = UK J) +UIND 
VKJ) = V1IJ)+VIN0 
W1IJ) = W1 ( J ) +WI ND 

221 CONTINUE 

PARI 6) = VEL JIT 
PAR I 5 ) = Fi 
PARI 9 ) = DIR1T 
DO 222 J*1,KSEG 

CALL VELOC 1 1 ,N3 , 21T.X1T, 0XDZIT,U J IT, D1T, UUE1T, XJ IT, YJ 1T,Z J IT, 

1 D JET IT ,CF1T,PAR,XP1 (J),YP1IJ),ZP1IJ ) ,UINO, VINO, WINO, SXZ1T ) 
UKJ) = U1IJJ-UIN0 
V1IJ) « VI I J ) -VINO 
WltJ) = W 1 I J ) —WINO 

222 CONTINUE 
250 CONTINUE 

RETURN 

ENO 


SUBROUTINE BITEST tI,TNEG,KS) 


TESTS FOR BLOCKAGE AND INTERSECTION, CALLED AS PART OF INTEGRATION 
LOOP 


DIMENSION Xl(ll,10),ZHll,l0),UJHll,10),01tll,10),DX0ZK 11,10) 
DIMENSION X21 11,10) ,Z2 1 11, 10) ,UJ2 1 11, 10 ), 02 1 11, 10 ) , 0X0Z21 1 1 , 10) 
DIMENSION X3I 100) .Z3I100) ,UJ3 ( 100 ) ,031 100 ) , DXDZ3 l 100 ) 

DIMENSION XBAS1I100) ,YBAS1(100) .ZBASIIIOO) 

DIMENSION XBAS2I 100) , YBAS2 (100 ) , ZBAS2I 100 ) 

DIMENSION XBAS3 1 100 ) ,YBAS3 1 100) ,ZBAS31 100 ) 

DIMENSION CF1(3,3, 10) ,CF2 1 3,3,10) ,CF1T1 3, 3 ) ,CF2T1 3, 3 ) ,CF3( 3,3) 
DIMENSION UUE1 1 11,10) ,UUE2 (11,10)«UUE1T 1100) ,UUE2T ( 100 ) ,UUE31 100) 
DIMENSION PARI15) 

DIMENSION XJ1(10),YJ1(10) , ZJ1 ( 10 ) , DJET1 ( 10 ) , VELJ 1(10) 

DIMENSION XJ2 ( 10) , Y J2 ( 10) , ZJ2 ( 10 ) , 0JET2 (10 ) , VEL J2 (10 ) 

DIMENSION ALFQl(lO) »BETQ1 ( 10) ,GETQ1( 10) , ALFQ2 ( 10 ) * BETQ2 ( 10 ) , 

1 GETQ2( 10) 

DIMENSION 0IR1(10),DIR2(10),ZS01(10),ZS02(10) , UFACT 1 ( 10), 
l UFACT2( 10) 
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COMMON /BLKl /CF 1 ,CF2 ,CFIT,CF2T ,CF3,UUE1, UUE2, UUE1T , UUE2T, UUE3,PAR 
C0MM0N/BLK2/X1 * Z1 »UJl ,01 , DXDZ1 , X2, Z2,UJ2, 02, DXDZ2 
C0MM0N/BLK4/X3»Z3»UJ3,D3,DXDZ3 

C0MM0N/BLK5/XBASI v YBAS1 * ZBAS1 , X8AS2, Y8AS2, ZBAS2, XBAS3, YBAS3, ZB A S3 
C0MM0N/BLK7/ALFQ, BETQ*GETQ*F1 ,F2»F3* VKONST 
C0MM0N/BLK8/ALFQ1 ,BETQ1 ,GETQ1 , ALFQ2,BETQ2, GETQ2 
C0MM0N/8LK9/MULT * IHOLD1 , KOUNT 1, IONE, I TWO* ITHR,N i*N2, N3, IF I XI 
C0MM0N/BLK10/XJ1,YJ1,ZJ1,DJET1,VELJ1,XJ2,YJ2,ZJ2,DJET2,VELJ2 
C0MM0N/BLK12/XJ3,YJ3,ZJ3,DJET3,VELJ3 
COMMON/BLK13/G,G2,G3,STEPI,STEPI2,STEPI3 
C0MM0N/BLK14/V2XI »V2Y1, V2Z1»V2X2»V2Y2* V2Z2 

C0MM0N/BLK15/DI Rl »0IR2»0IR1T ♦ DI R2T ,DR3, ZSOly ZS02» UFACT1* UFACT2 
C0MM0N/BLK17/GS,A,B.C*IS1.IS2,NPS 
C 

DE * .OOOl*OJETUl) 

IF (MULT-2) 21,200,200 

200 IF ( I H0LD1-1 ) 201,202,21 

201 IF ( TNEG) 203,203,204 

203 CALL XPROO ( V2X1 , V2Y1 , V2Z1 ,ALFQ1 ( KS ) , BETQ1 (KS ) , GETQ1 (KS ) , XT 1, 

1 YT1 , ZT1 ) 

CALL XPROO (XT1,YT1,ZT1 , ALFQ1 (KS) , BETQ1 (KS ) ,GETQ1 (KS),CFNX, 
i CFNY,CFNZ) 

CALL PLANE (CFNX,CFNY,CFNZ,XBAS1 (I) ,YBAS1 ( I),ZBASl( I ),V2X2,V2Y2, 

1 V2Z2 ,XJ2 (KS ) ,YJ2(KS) , Z J2 ( KS ) ,X INT , Y INT, Z I NT ) 

IF ( Y INT— YJ2 ( KS )— DE ) 205,205,22 

204 UUE2( I ,KS) * 1. *UFACT2 ( KS > 

CALL XPROO ( V2X2,V2Y2,V2Z2,ALFQ2(KS) , BETQ2 (KS ) , GETQ2 ( KS ) , XT2, 

1 YT2,ZT2> 

CALL XPROO (XT2,YT2,ZT2,ALFQ2(KS), BETQ2 (KS ) , GETQ2 ( KS ) *CFNX* 
l CFNY ,CFNZ) 

CALL PLANE ( CFNX,CFNY,CFNZ,XB AS2 ( I } , YBAS2 ( I ) , ZBAS2( I ) , V2X1 , V2Y1, 

1 V2Z1 ,XJ1 (KS) , Y Jl ( KS ) , Z J1 ( KS ) ,XINT,YINT,ZINT ) 

IF (YINT-YJK KS )— DE ) 205,205,22 

205 I HOLD 1 = 1 

202 IF (TNEG) 206,206,207 

206 I TWO = I-K0UNT1 
GO TO 208 

207 I ONE * I-KOUNTl 

208 IT1 = I ONE 
IT2 = I TWO 
N1 = IT1+1 
N2 = IT2+1 

IF (KS-U 210,210,211 

210 I SI = I ONE 
IS2 = I TWO 
GO TO 212 

211 IS1 = I- ( KS-1 ) *NPS 
IS2 = ISi 

212 CONTINUE 

CALL COMP (V2X1,V2Y1,V2ZI,V2X2,V2Y2*V2Z2,XBAS1( IT1),YBAS1( IT1), 

1 ZBAS1( IT1! ,XBAS2( IT2),YBAS2( IT2 ) , ZBAS2 ( I T2 ) , Z1 1 I SI ,KS ) ♦ 

2 Z2( IS2,KS),D1( ISi,KS) ,DJET1 (KS),02( IS2,KS ) ♦ 0JET2 ( KS ) ,VELJ1(KS) , 

3 VELJ2(KS),DX0Zi( ISI ,KS ) ,KS,UUE2( I S2 , KS ) , Al, A2, 0R3, INT) 

IF (INT) 21,21,209 

209 IH0LD1 = 2 
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Ml = ITl 
N2 = IT2 
IFIX1 = I 

CALL BALANC < XBASU ITl ) , YBAS1 ( ITl ) , ZBASK ITl ) , XBAS2 ( IT2 ) , 

1 YBAS2 (IT2),ZBAS2(IT2)«UJ1(IS1,KS),UJ2( IS2,KS ) * VELJ l (KS ) • 

2 VELJ2(KS),A1,A2*V2X1,V2Y1,V2Z1*V2X2*V2Y2*V2Z2*DR3»XJ3»YJ3*ZJ3* 

3 DJET3* V2X3, V2Y3* V2Z3 *VEL J3) 

CALL CFCAL ( ALFQ*BETQ*GETQ,V2X3» V2Y3* V2Z3,CF3 ) 

CALL ROTATE ( V2X3,V2Y3,V2Z3,CF3,VXT,VYT,VZT,0) 

U J3 ( 1 ) = l. 

03(1) = l. 

X 3 ( 1 ) = 0. 

Z3 ( l ) * 0, 

DXDZ3 ( 1 ) » VXT/VZT 
XBAS311) * XJ3 
YBAS3 ( l ) * YJ3 
ZBAS311) * ZJ3 
D = ATAN( VXT/VZT ) 

IF ( VXT) 901,902*902 

901 F3 = .3*CGS(0) 

GO TO 903 

902 F3 = .3/C0S(D) 

903 CONTINUE 

G3 = GS+DJETl ( 1 ) /DJET3 
STEPI3 = .2*G3 
GO TO 21 

22 KOUNTl = KOUNTl+i 

151 = I 

152 = I 
21 CONTINUE 

RETURN 

ENO 


SUBROUTINE INTEG ( I ,TNEG*KS) 

INTEGRATION OF THE EQUATIONS OF MOTION FOR THE JET PATH 

EXTERNAL DERIV 
C 

DIMENSION XI ( 11,10),Z1(11*10) *UJ1 ( 11, 10 ) , 01 ( 11, 10 > * DXDZ 1< 11,10) 

01 MENS I ON X2(11,10),Z2(11,10) ,UJ2( 1 1 , 10 ) , 02 ( 1 1, 10 ) , DXDZ2( 1 1 ♦ 10 ) 
DIMENSION X3( 100) ,Z3(100) ,UJ3 ( 100 ) ,03 ( 100 ) , DXDZ3( 100) 

DIMENSION XBASU 100) ,YBASl (100) ,ZBAS1 ( 100 ) 

DIMENSION XBAS2 ( 100 ) , YBAS2 ( 100) , ZB AS 2 ( 100) 

DIMENSION XBAS3( 100) ,YBAS3(100) , ZB AS 3(1 00) 

DIMENSION CF1(3,3,10) ,CF2 ( 3,3 , 10 ) , CFIT ( 3, 3 ) ,CF2T ( 3, 3 ) ,CF3( 3,3) 
DIMENSION UUEK 11,10) ,UUE2( 11,10 ),UUE1T( 100), UUE2TI 100 ) ,UUE3( 100) 
DIMENSION SDXDZKll.lO) ,SDXDZ2(11, 10 ) , SXZ IT ( 100 ) , SXZ2TJ 100), 

1 SDXDZ3( 100) 

DIMENSION PARI 15) 

DIMENSION XJU10),YJ1(10> , ZJl ( 10) , DJET1 ( 10 ) , VELJ l ( 10 ) 

DIMENSION XJ2(10) ,YJ2(10) , ZJ2 ( 10 ) , DJET2 ( 10 ) , VEL J 2 ( 10) 
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DIMENSION ALFQl(iO) ,BETQi(10) ,GETQl 1 10) , ALFQ2( 10) ,BETQ2 ( 10) , 

1 GETQ21 10) 

DIMENSION DIR1(10)«DIR2(10),ZS01(10),ZS02(10 ) »UFACTI ( 10) * 
l UFACT2( 101 

C0MM0N/8LKI/CF1,CF2,CF1T,CF2T,CF3,UUEI,UUE2,UUE1T,UUE2T,UUE3,PAR 
COMMON /8LK2/X 1 , Z1 *UJ1 *D1 * DXDZ1 «X2« Z2,UJ2, D2, DXDZ2 


C0MM0N/BLK4/X3* Z3*U J3 *D3 * DXDZ3 

COMMON /BLK5/XBAS1 , YBAS1 , ZBASI , XBAS2, YBAS2, ZBAS2, XBAS3, YBAS3, ZBAS3 
C0MM0N/BLK7/ALFQ, BETQ,GETQ,Fl,F2,F3,VK0NST 
C0MM0N/BLK8/ALFQ1 tBETQI ,GETQ1 « ALFQ2, BETQ2«GETQ2 
COMMON/BLK9/MULT,IHOL01,KOUNT1,IONE, ITWO, ITHR»Ni»N2,N3* IFIXl 
COMMON/BLKIO/X J1 ,YJ1,ZJI« DJET 1 , VEL J1 , XJ2* YJ2, ZJ2, DJET2, VEL J 2 
C0MM0N/BLK12/XJ3, YJ3»ZJ3»DJET3»VELJ3 
C0MM0N/8LK13/G*G2 *G3» STEPI »STEPI2*STEPI3 
C0MM0N/8LK14/V2X1 * V2Y1 * V2Z1 * V2X2, V2Y2* V2Z2 

C0MM0N/BLK15/DI Ri *DIR2 1 DIR1T « DIR2T , DR3, ZS01, ZS02* UFACT1, UFACT2 
C0MM0N/BLKI6/SDXDZ1 ,SDXDZ2,SDXDZ3, SXZ1T*SXZ2T 
C0MM0N/BLK17/GS ,A,B ,C, I SI * IS2«NPS 

DIMENSION FIN(4),F0UT(4) 

IF (MULT-21 24,51,51 
51 IF ( IH0L01— 2 J 25,30,30 
25 IF (TNEGJ 24,24,27 

27 IF (IH0LD11 28,28,24 

24 UUE1 ( IS1 «KS) = l.*UFACTl(KS) 

PARI 6 ) * VEL J 1 ( KS) *UUE1 ( IS1,KS) 

PARI 5 ) = FI 
P AR( 9 ) = OIRl(KS) 

P AR( 15) * 1. 

Z l( ISl+l«KS) = Zlt IS1«KS)+G 
F IN( 1 1 = UJ1 ( I SI ,KS ) 

F IN( 2 ) * Dl ( I SI ,KS ) 

F IN( 3 1 = XI ( I SI ,KS) 

F IN( 4 1 = DXDZ 1 ( I SI ,KS) 

CALL ADAMS ( 4,Z1 ( I SI ,KS ) , Z1 ( I Sl + 1 ,KS ) , STEPI , G, 999, 1 .0E-04, 1 .0E-05, 
i 0,FIN*F0UT , PAR,DERIV 1 
UJ1( ISl+l,KS) = FOUTll) 

Di ( IS1 + 1 ,KS) = FOUT ( 2 ) 

XI ( I Sl+1 »KS) = FOUT (3) 

DX0Z1 ( ISl+l*KS) = FOUT ( 4 ) 

SDXDZK ISl+l*KS) * PARI 10) 

CALL 0UTPT1 (Xl(ISL+l,KS)«Zl(ISl+l«KS)«DXDZl(ISl+l,KS)«CFl,KS, 

1 0JET1 ( KS) ,XJ1(KS),YJ1(KS),ZJ1(KS) , XBASl ( IONE+1 ), YBASlt IONE+1), 

2 ZBASU IONE+1), V2X1,V2Y1,V2Z1) 

IF (MULT-2) 50,41,41 

41 IF (IH0L01) 50,50,28 

28 PAR( 6 ) = VEL J2( KS) *UUE2 ( IS2,KS) 

PAR( 5 ) = F2 

PARC9) = DIR2 (KS) 

PAR( 15) = 2. 

Z2( I S2+1 *KS ) = Z2( IS2,KS)+G2 
FIN(l) = UJ2 ( I S2 ,KS ) 
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FIN(2) = D2( IS2»KS) 

FIN<3> = X2( I S2 «KS) 

FIN(4) = 0XDZ2! IS2.KS) 

CAU ADAMS ! 4 *Z2 ( I S2 .KS >. 12 ( I S2 + I.KS ), STEPI2, G2.999, 1 .0E-04, 

1 1.0E-05.0, FIN, FOUT, PAR, DERIV) 

UJ2! IS2+I.KS) * FOUT(l) 

D2 ( I S2+1 »KS ) * FOUT ( 2 ) 

X2 ( IS2+1 »KS ) * FOUT l 3 ) 

0XDZ2( IS2+l«KS) * FOUT (4) 

SDXOZ2 ( I S2+1 *KS ) * PAR(IO) 

CALL 0UTPT1(X2( IS2+1 ,KS ) , Z2 (IS2+1 ,KS ) , DXDZ2 ! I S2+1 , KS ) ,CF2, KS, 

1 DJET2! KS ) » X J2 ( KS ) , Y J2! KS ) ,Z J2 ( KS ) * XBAS2 ! ITWO+1 ) * YBAS2 ( ITWO+1), 

2 ZBAS2!ITWO+l),V2X2»V2Y2»V2Z2) 

GO TO 50 

30 I THR = I-IFIXi+l 
UUE3( I THR) * 1. 

PAR ( 6 ) = VELJ3 
PAR ( 5 ) = F3 
PAR ( 9 ) = 0R3 
PAR (15) = 3. 

Z3IITHR+1) * Z 3 ( ITHR ) +G3 
F IN ( 1 ) = UJ3 ( ITHR) 

F IN( 2 ) = D3! ITHR) 

F INI 3 ) = X3IITHR) 

FINU) = DX0Z3( ITHR) 

CALL ADAMS! 4, Z3< ITHR) ,Z3( ITHR+1 ) ,STEPI3,G3,999, 1.0E-04, 

1 l«0E-05 ,0»F I N, FOUT , PAR, DERI V ) 

UJ3UTHR + 1) = FOUT ( l ) 

D3< ITHR+l ) * FOUT { 2 ) 

X3( ITHR+1 } = FOUT ( 3 ) 

DXDZ3 ( I THR+1 ) = FOUTU) 

SDXDZ3! I THR+ l ) = PAR(IO) 

CALL OUTPT ( X3< ITHR+1 ) * Z3 ( ITHR+1 ) , DXDZ3 ( ITHR+1 ),CF3, DJET3, XJ3, YJ3, 
1 ZJ3.XBAS3! ITHR+1 ),YBAS3! ITHR+1), ZBAS3! ITHR+i ), V2X3, V2Y3, V2Z 3 > 

50 CONTINUE 
RETURN 
END 


SUBROUTINE DERI V (Z,FN,FPR,PAR) 

COMPUTES DERIVATIVES FOR ADAMS PREDICTOR/CORRECTOR METHOD 

DIMENSION FN( 1 ) «FPR ( 1 ) f PAR ( 1 ) 

C 

El = PARll) 

E2 = PAR ( 2) 

E 3 = PAR ( 3) 

F * PAR ( 5) 

VEL J=PAR ( 6) 

PI » PAR! 7 ) 

Cl = PAR! 8) 

DR = PAR ! 9) 

UJ = FN ( 1 ) 
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D * FNC2) 

DXDZ*FN(4) 

JET * PAR ( 15) 

IF (JET-2) 1,2,5 

1 A = PARC 11) 

B = PAR (12) 

GO TO 3 

2 A = PARI 13) 

B = PARC 14) 

3 CONTINUE 

IF (A) 5,5,4 

4 IF CZ-A/B) 6,6,7 

6 AN1 = (A*A-2.*A*B*Z+B*B*Z*Z)/(D*D) 

AN2 = PI*UJ*CA*B-B*B*Z)/2. 

AN3 = 1.-4. *ANl 
GO TO 10 

7 IF CJET-l) 8,8,9 

8 PAR(ll) = 0. 

GO TO 5 

9 PARC 13) = 0. 

5 AN1 = 0. 

AN2 » 0. 

AN3 = 1. 

10 CONTINUE 

COST a 1. /SORT C l.+DXDZ*DXDZ) 

SINT a SIGN! 1 • ,OXDZ ) *SQRT C l.-COST*COST ) 

E = E2/C l«*E3*C0ST/( VELJ*UJ) ) 

IF (VELJ*UJ-SINT) 11,12,12 

11 E * o. 

12 ZSO * C 1.— DR) *VELJ*F/.75 
ZP = Z+ZSO 

IF C ZP— VELJ*F ) 47,60,60 

47 IF CZP-10.) 40,60,60 
40 IF (ZP-.6*VELJ*F) 42,43,43 

42 E a E*. 1/.32 
GO TO 60 

43 IF CZP-.8*VELJ*F) 44,45,45 

44 E = E*. 12/. 32 
GO TO 60 

45 E * E*.21/.32 
60 ZOVM a ZP/VELJ 

IF CZOVH-F) 22,23,23 

22 VAR = SQRT((l.-Ml.-.75*Z0VM/F)**2)/2.) 

XT » 1.-.754Z0VM/F 
XT « 1. /XT 

CO = C-XT*XT+6.6*XT+.4>/6. 

VAR1 a E1*C0ST+E*(VELJ*UJ-SINT)*P1*VAR 
VAR2 * VEL J*VEL J*C0ST 

VAR3 = . 25*PI*( l.-.75*Z0VM/F— AN1 )*UJ*D 

VARA = .25*PI*C1.-.75*Z0VM/F)*UJ*0 

DUJ = C VAR1*SINT/VAR2— VAR1*U J/ ( VELJ*C0ST ) )/VAR3 

DO = CVARl*0/CVELJ*C0ST)+3.*PI*0*0*UJ/C16.*F*VELJ)-VAR3*O*0UJ/ 
l UJ-AN2 ) /C 2. *VARA) 

VAR4 = C E1*.5*C0 ) *C0ST+E*C VEL J*UJ— SINT ) *Pl*VAR 
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DDXDZ= VAR4/I VAR2*C0ST*VAR3*UJ ) 

GO TO 15 

23 VAR1 = El*CO$T*E*<VELJ*UJ-SINT)*Cl 
CD = 1.8 

DUJ = 16.*VAR1*( SI NT/I VEL J*VEL J*COST)-UJ/ 1 VELJ*COST ) )/IPI*D*UJ ) 
DUJ » DUJ/AN3 

DO = 8.*( VARl/<VELJ*C0ST)-AN3*PI*D*DUJ/16.-AN2/0)/(PI*UJ) 

VAR4 = I E1+.5*CD) *COST+E*(VELJ*UJ-SINT)*Cl 
DDXDZ® 16.*VAR4/ 1 PI *VEL J*VEL J*D*UJ*UJ*COST*COST ) 

DDXDZ = DDXDZ/AN3 
15 CONTINUE 

PARI 10 } = DDXOZ 
FPR(l) = DUJ 
FPR( 2 ) = OD 
FPR ( 3 ) = OXDZ 
FPRI4 ) = DDXDZ 
RETURN 
END 


SUBROUTINE VELOC ( N1 * N2 « Z * X , DXDZ,U J » D,UUE* XJ, YJ , Z J, DJET,CF, PAR * 
1 X0*Y0*Z0*UIF*VIF*WIF,D2XDZ2) 

EVALUATES INDUCED VELOCITIES AT ONE CONTROL POINT (XO.YO.ZO IN 
FIXED COORDINATE SYSTEM) FOR A GIVEN JET 

C0MM0N/8LK19/DI ARAT »DREF 
C 

DIMENSION ZI 1) ,Xt 1) ,DXDZI 1 ) »U J 1 1 ) * Dl 1 ) , UUE I 1 ),PAR( 1 ) 

DIMENSION CF I 3* 3 ) 

DIMENSION D2XDZ2 ( 1 ) 

C 

E2 = PARI 2) 

E3 = PARI 3) 

F = PARIS) 

VEL J=PARI 6) 

PI = PARI7) 

Cl * PARI 8) 

DR = PARI9I 
N = N2-N1+1 
IF I N/2-I N+l ) /2 ) 1*2,2 

1 M = (N-l)/2 
GO TO 3 

2 M * I N-2 ) /2 

3 XPT » I XO-XJ ) /D JET 
YPT = <Y0-YJ)/0JET 
ZPT * ( ZO-Z J ) /DJET 

CALL ROTATE I XPT ,YPT*ZPT»CF*A*6*C*0) 

UI = 0. 

VI = 0. 

WI = 0. 

Ml = M+l 
00 21 K=N1*M1 
El = PARI 1 ) 
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IF (K-M) 11,11,10 

10 IF (N/2-(N*lJ/2) 22,12,12 
12 I * 2*K-l 

ZINCR = Zd+lJ-Zd) 

GO TO 14 

11 I * 2*K 

ZINCR = ZU + 1J-ZU-1) 

14 COST = l./SQRT(l.«-DXDZ<I>*OXDZUn 

SINT = SIGN(1.,DXDZ<I))*SQRT(1.-C0ST*C0ST> 

S IE = -( (Z(I >-C)*C0ST+(X( I )-A)*$INT) 

ETA = B 

ZETA= ( Z ( I >-C ) ♦SINT- 1 X C I )-A)*COST 
01 = .5*0(1) 

DCUB1 = SIE*SIE+ETA*ETA+ZETA*ZETA 
D0UB2 = SORT ( D0UB1 ) 

UBLOCK = .5*D1*DI*ZI NCR*COST* < 1 .-3 »*ZET A*ZET A/D0UB1 ) / ( DOUB 1*D0UB2 ) 
1 -SINT*1.5*SIE*ZETA*Dl*Dl*ZlNCR/( D0UB1*D0UB 1*D0UB2 ) 

VBLOCK = -1.5*ZETA*ETA*D1*D1*ZINCR/(DOUB1*OOUBI*DOUB2) 

WBLOCK = -.5*D1*01*ZINCR*SINT*(1.-3.*ZETA*ZET A/ DOUB 1 ) / ( DOUB l* 
l D0UB2)— C0ST*l.5*SIE*ZETA*Di*Dl*ZINCR/( DOUB 1*D0UB 1*D0UB2 ) 

VELJE = VEL J*UUE ( I ) 

CURV = 02X0Z2 ( I ) / ( ( 1. +DXDZ ( I ) *DXDZ < l ))**1 .5 ) 

CURV = 3.*CURV*0REF/DJET 

El = El-CURV/COST 

E * E2/(1.+E3*C0ST/(VELJE*UJ(I))) 

IF (VELJE *11 J( I ) -SI NT) 51,52,52 

51 E = 0. 

52 ZSO = (l.-DR)*VELJE*F/.75 
ZP = Z(I)+ZSO 

IF (ZP-VELJE*F) 47,60,60 
47 IF (ZP-10.) 40,60,60 
40 IF (ZP-.6*VELJE*F> 42,43,43 

42 E = E*. 1/.32 
GO TO 60 

43 IF l ZP— . 8*VEL JE*F ) 44,45,45 

44 E = E*. 12/. 32 
GO TO 60 

45 E = E*.21/.32 
60 ZOVM = Z P/VELJE 

IF ( ZOVM— F ) 31,32,32 

31 VARB = < l.-.375*Z0VM/F) 

VAR = SQRT( ( 1«+(1.-.75*Z0VM/F) **2 ) / 2. ) 

HT3 = .25*ZINCR*(E1+E*PI*VAR*(VELJE*UJ(I)-SINT J/COST) 

GO TO 33 

32 VARB = .625 

HT3 = .25*ZINCR*(E1+E*t VELJE*UJ( I )-SINT)*Cl/COST) 

33 UBLOCK * UBLOCK*VARB 
VBLOCK = VBLOCK*VARB 
WBLOCK = WBLOCK+VARB 

Z1 * (C-Zm)*(C-ZU)) + (A-X(I))*(A-X<IJ) 

Z2 = SORT 1(B-D1)*(B-D1)*Z1) 

Z 3 = SORT ((B+D1)*(B+D1)+Z1) 

USINK = -HT3* ( X ( I ) -A) *1 (B-Dl )/(Zl*Z2)-(B*Dl)/(Zl*Z3) )/PI 
VS INK = -HT3*< 1./Z2-1./Z3J/PI 
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WSINK * -HT3*(ZU )-C> *1 (B-Dl ) / < Zl*Z2 l-< B+Dl >/ ( Z14Z3 M /P I 
IF ( UUE( I )-l. ) 6,5,6 
6 FACT = l./UUEd ) 

UBLOCK a UBLOCK*FACT 
VBLOCK = VBLOCK*FACT 
UBLOCK » WBLOCK*FACT 
USINK = USINK*FACT 
VSINK = VSI NK*FACT 
WSINK = HSINK*FACT 
5 U! * U I +U SINK ♦UBLOCK 
VI = VI ♦VSI NK+ VBLOCK 

21 WI = WI+WSINK+W BLOCK 

22 CALL ROTATE < UI F , VI F , WI F,CF,UI , VI , W l , 1 ) 

691 FORMAT (6F12.5) 

RETURN 

END 


SUBROUTINE COMP ( VXl , VY1 , VZI, VX2 . VY2. VZ2.X 1 . VI . 7 1 . X?. V7. 7 7. 7 1 I t 7?i t 
1 01 ,OJl ,02,0J2 ,VI,V2,SL1,KS,UUEFF, Al, A2, DRAT » IND ) 

COMPUTES U/UEFFECTIVE AND TESTS FOR INTERSECTION OF CENTERLINES 

DIMENSION DIRK 10 > ,DIR2 ( 10 ) , ZSOK 10) , ZS02 ( 10 ) ,UFACT1 ( 10 ) , 

1 UFACT2 (109 

C0MM0N/BLK7/ALFQ,BETQ,GETQ,F1 ,F2» F3, VKONST 

C0MM0N/BLK15/DI R1 ,DIR2,DIR1T, DIR2T , DR3, ZS01,ZS02*UFACT1»UFACT2 
C0MM0N/BLK19/DI ARAT,DREF 

IND a 0 
PI = 3.1416 

CALL XPROO ( VXl ,VY1 , VZI « ALFQ, BETQ, GET Q, CFNX, CFNY* CFNZ ) 

CALL XPROO ( VX2,VY2,VZ2, ALFQ, BETQ, GETQ, XT2, YT2, ZT2 ) 

CALL PLANE < CFNX , CFNY, CFNZ , XI , Yl , Z1 , XT2, YT2, ZT2, X2, Y2, Z2, X I , YI , Z I ) 
DIST = SQRT( < XI-X2) **2+I YI-Y2 )**2+ (ZI— Z2)**2 ) 

COMPUTE U/UEFFECTIVE 

R = D1*DJ1*.5-DIST 

FACT = (i.0+R/(D2*DJ2*.5) >*.5 

IF C FACT-1. ) 10,10,11 

11 UUEFF = VKONST 
GO TO 15 

10 IF ( FACT ) 13,13,12 

13 UUEFF = 1. 

GO TO 15 

12 TEST1 = Dl*DJl*. 5+DI ST 
TEST2 * D2*0J2*.5 

IF (TEST1-TEST2) 16,14,14 
16 FACT = D1*DJ1/(D2*DJ2) 

14 UEFU = 1.4-fl. /VKONST-l. )*FACT 
UUEFF = l./UEFU 

15 CONTINUE 
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UUEFF = UUEFF*UFACT2 ( KS ) 

TEST FOR INTERSECTION OF CENTERLINES 

COST = l • /SORT (1.+SL1*SL1) 

SUMO = 0 J 1*01 *• 5 
IF (DIST-SUMD) 22,99,99 

22 OISTN = SQRT ( < Xl-XI ) **2 + ( Yl-Yl ) **2*( Zl-ZI )**2 ) 
IF ( 0IR1 ( KSJ-.2501 ) 25,25,51 

51 ZOVM = tZS01(KS)*ZlL)/(Vl*UFACTl(KSn 
IF { ZOVM-Fl ) 24,24,25 

24 FACT1 * l.-.75*Z0VM/Fl 
GO TO 26 

25 FACT1 = .25 

26 IF <DIR2(KS)-.2501) 28,28,52 

52 ZOVM = ( ZS02 < KS ) +Z2L ) / ( V2*UUEFF ) 

IF ( Z0VM-F2 ) 27,27,28 

27 FACT2 = 1 7 5* ZOVM /F2 
GO TO 29 

28 FACT2 = .25 

29 SUMO = DJ1*D1*FACT1*C0ST*.5 
IF (DISTN-SUMO) 30,30,40 

30 IND * 1 
GO TO 45 

40 IF ( X2-X1 ) 30,30,99 

45 A 1 = PI*FACTi*01*Dl*DJl*DJl*.25 
A2 = PI*FACT2*D2*D2*DJ2*DJ2* .25 
DRAT = OIARAT 

99 CONTINUE 
RETURN 
END 


SUBROUTINE VEL1 < MULT , ALFA, VK1 ) 

COMPUTES EFFECTIVE VELOCITY RATIO FOR DOWNSTREAM JET AT EXIT 

DIMENSION XJK10) ,YJ1(10> , ZJ1 < 10 ) , DJET1 ( 10 ) , VEL J It 1 0 ) 

DIMENSION XJ2(10) ,YJ2(10) ,ZJ2 ( 10 ) , DJET2 ( 10 ) , VEL J2 ( 10) 

C0MM0N/BLK7/ALFQ,BETQ,GETQ,F1 ,F2,F3,VK0NST 

COMMON /BLK1 0/XJl»YJl,ZJl»DJETl,VELJl»XJ2*YJ2»ZJ2,DJET2»VELJ2 
COMMON /BLK14/V2 XI ,V2Y1,V2Z1,V2X2,V2Y2,V2Z2 

VELJl(l) » l./VELJKl) 

IF (MULT-21 15,1,1 
1 VELJ2 ( 1 ) = 1 . /VEL J2 ( 1 ) 

DOTP = (XJ2( 1)-XJ1 ( 1 ) )*ALFQ+(YJ2U)-YJ1(1) )*BETQ + (ZJ2( 1 J-ZJll 1 ) ) 
l *GETQ 

DEN = SQRT((XJ2(l)-XJl(l))**2 + (YJ2(i)-YJl(l))**2-MZJ2(l)-ZJl(l) ) 
l **2) 

DOTP = DOTP/OEN 
IF ( ABS(D0TP)-.02) 10,10,11 
10 VK1 = 1. 

GO TO 15 
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11 CONTINUE 

A * ALFA*. 0174533 
ALF = COS( A) 

BET * SINCA) 

GET * 0. 

CALL XPROO IV2X1,V2Y1,V2Z1,ALF,BET,GET,XT1,YT1,ZT1) 

CALL XPROO < XT1 » YT1 »ZTl » ALF,BET ,GET,CFNX, CFNY»CFNZ ) 

CALL PLANE ( CFNX,CFNY,CFNZ ,X J1 ( 1 ) , Y J1 ( 1 ) . ZJ1< 1 ) , V2X2, V2Y2, V2Z 2, 
1 XJ2( 1 ) . Y J2 ( 1 ) * Z J2 ( I ) *XI,YI,ZI) 

S = SORT (IXJl(l)-XI) **2+ ( Y J1 ( l)-YI)**2+(ZJl(l)-ZI)**2)/DJETi(l) 
VK1 = I S+.75 ) /< S— 1 • ) 

OOTP = V2X1*ALFQ*V2Y1*BETQ+V2Z1*GETQ 
AIN * ACOSIDOTP) 

IF IDOTP) 4,4,5 

4 AIN = AIN-3.14159/2. 

GO TO 6 

5 AIN = 3. 14159/2. -AIN 

6 CONTINUE 

SIN2 * SIN(AIN)*SIN(AIN) 

C0S2 = COS( AIN) *COS ( AIN ) 

C = 1./VK1 

VK1 * l./SQRT(SIN2*C*C*C0S2) 

15 CONTINUE 
RETURN 
END 


SUBROUTINE BALANC ( XI , Y1 , Z1 ,X2 , Y2, Z2,UJ 1, UJ2, VI, V2, Al, A2, VX1 , VY1 , 

1 VZi,VX2,VY2,VZ2, FACT 1,X3«Y3«Z3*DJ3,VX3,VY3,VZ3, 

2 VELJ3) 

ESTABLISHES INITIAL CONDITIONS FOR NEW JET FROM MOMENTUM BALANCE 

PI = 3.1416 
X3 = (Xl+X2)*.5 
Y3 = ( Y1+Y2) *.5 
Z3 = (Z1+Z2) *.5 
XMi = UJI*V1*A1 
XM2 = U J2*V2*A2 
DEN * XM1+XM2 

UJX * ( XM1*UJI*V1*VX1+XM2*UJ2*V2*VX2)/DEN 
UJY * ( XM1*UJI*V1*VYH-XM2*UJ2*V2*VY2)/DEN 
UJZ * { XM1*U J1*V1*VZ1*XM2*UJ2*V2*V Z2 )/DEN 
VELJ3 = SORT <UJX*UJX+UJY*UJY+UJZ*UJZ) 

VX3 * UJX/VELJ3 
VY3 * UJY/VELJ3 
VZ3 * UJZ/VELJ3 
A3 * DEN/VEL J3 

DJ3 » SORT <4.*A3/IPI*FACT1)) 

RETURN 

END 
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SUBROUTINE FIX ( UA,VA,WA,UB,VB,WB,NO) 

LIMITS MUTUALLY INDUCED VELOCITIES TO A MAXIMUM VALUE 
DIMENSION UA(l) ,VAU),WA< I ) ,UB( l ) , VB( 1 ) * WB ( I ) 

DEL = .6 



DO 5 I -1 ,N0 

IF ( ABS < UAI I ) )-DEL) 

1,1,3 

1 

IF (ABS(VA(D)-DEL) 

2,2,3 

2 

IF (ABS< WA( I ) J-DEL) 

5,5,3 

3 

UAU) = UA(I-l) 


5 

VA(I) = V A ( 1 — 13 
WA(I) = WAU-1) 
CONTINUE 



DO 15 1=1, NO 
IF (ABS(UB(I) )-DEL) 

11,11,13 

11 

IF ( ABS( VB( I ) ) -DEL ) 

12,12,13 

12 

IF ( ABS ( WB( I ) )-DEL ) 

15,15,13 

13 

UB(I) = UB(I-l) 


15 

VB ( I ) = VB(I-l) 
WB ( I ) = WB(I-l) 
CONTINUE 



RETURN 

ENO 



SUBROUTINE MUINT ( ALF , BET ,GET , U, V, W, ALFM, BETM.GETM, UUIND ) 

COMPUTES DIRECTION COSINES OF MODIFIED FREESTREAM AND THE 
MODIFIED FREESTREAM/FREESTREAM VELOCITY RATIO 

A = ALF+U 
B = BET+V 
G = GET+W 

D = SQRT ( A**2«-B**2+G**2 ) 

ALFM = A/D 
BETM = B/D 
GETM = G/D 
UUIND = l./D 
RETURN 
END 


SUBROUTINE PRTOUT ( IGEOM,XO,YO,ZO,U*V,W,CP*NK* ITR, ITF,OV) 

PRINTS OUT COMPUTED ANSWERS, INFORMATION INCLUDES JET CENTERLINE 
OATA AND INDUCED VELOCITIES AT CONTROL POINTS 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 


XIT<100),Z1T(100),UJ1TUOO),D1T(100),DXDZIT( 100) 
X2T(10O)»Z2T(I00)tUJ2T(lOO)tD2T(IOO)» DX0Z2T (100) 
X3< LOO) ,Z3( LOO) , U J3 ( 1 00 ) , 03 ( LOO ) , OXDZ 3 ( 100) 

XBASl ( 100 ) ,YBAS 1(100) ,ZBAS 1(100) 

XBAS2 ( 100 ) , YBAS2 (100),ZBAS2(100) 

XB A S3 (100) ,YBAS3(100),ZBAS3( 100) 
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C0MM0N/BLK3/X1T,Z1T,UJ1T»DIT,DXDZIT,X2T,Z2T,UJ2T,D2T,DXDZ2T 

COMMON/BLK4/X3,Z3,UJ3,03,DXOZ3 

C0MM0N/BLK5/XSAS1 * YBAS1 • ZBAS1 ,XBAS2, YBAS2 , ZBAS2, XBAS3, YBAS3,ZBAS3 
COMMON /BLK9/MULT « IHOLD1 *KOUNTl, I ONE* I TWO* ITHR,N1,N2,N3, IFIXl 
COMMON/BLK12/XJ3, YJ3»ZJ3*DJET3* VEL J3 
C 

dimension xom,Yom,zom,u(i),vm,wm,cp<i) 
c 

WRITE (6*640) 

IF (OV-.Ol) 15,15,10 
10 WRITE (6,601) 

601 FORMAT (1H0,//) 

IF (ITR) 5,5,6 

5 WRITE (6,606) 

606 FORMAT ( 1H0,45X ,29H*** INITIAL APPROXIMATION *** ) 

GO TO 15 

6 IF (ITR-ITF) 7,8,8 

7 WRITE (6,607) ITR 

607 FORMAT ( 1H0,47X ,20H*** ITERATION NUMBER, 1 2, IX, 3H*** ) 

GO TO 15 

8 WRITE (6,608) ITR 

608 FORMAT ( 1H0, 39X ,20H*** ITERATION NUMBER , 1 2, 21H, FINAL ITERATION ** 
1*) 

15 WRITE (6,601) 

IF (MULT-2) 1,2,2 

1 WRITE (6,602) 

602 FORMAT l IH0,46X ,27H** SINGLE JET CENTERLINE **) 

GO TO 20 

2 WRITE (6,603) 

603 FORMAT ( 1H0,43X,33H** CENTERLINES OF JETS l AND 2 **) 

IF ( IH0LD1-2 ) 20, A, A 

A WRITE (6,605) 

605 FORMAT ( 1H ,51X,17HAND COALESCED JET) 

20 CONTINUE 

WRITE (6,630) 

630 FORMAT ( 1H0,42X,35H***********************************// ) 

IF ( MULT.GE. 1 ) WRITE (6,610) 

IF (MULT.GE. 2) WRITE (6,611) 

610 FORMAT ( 1H0, 3X,6HXCCORO,3X,6HYCOORO, 3X, 6HZC00RD, 3X, 2HUJ ,4X, 3HDI A ) 

611 FORMAT ( 1H+, 42X«6HXCG0RD« 3X,6HYC00RD» 3X ,6HZC00RD, 3X» 2HUJ, AX, 3HDI A ) 
WRITE (6,612) 

612 FORMAT (1H0) 

IF (MULT-2) 30,40,40 
30 CONTINUE 

WRITE (6,616) (X8AS1(I)«Y8AS1(I) , ZBAS1 ( I ) ,UJ 1T( I ) , D1T ( I), 1*1, Nl) 
616 FORMAT ( 1H , 1X,F8.2 , IX, F8.2, IX, F8. 2, IX, F5 . 3, IX, F5.2 ) 

GO TO 90 

40 IF (N1-N2) 41,42,42 

41 IP1 = Nl 
IP2 * N2 
GO TO 43 

42 I PI * N2 
IP2 = Nl 

43 CONTINUE 
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DO 47 1=1, IP1 

47 WRITE (6,613) XBAS1 ( I ) • YBAS1 ( I), ZB AS 1(1 ),UJ1T( I), D1T (I),XBAS2(I), 
1 YBAS2< I ) ,ZBAS2( I ) «UJ2T ( I ) ,D2T(I) 

613 FORMAT ( 1H , IX, F8.2 , IX, F8 .2, IX, F8. 2, IX, F5 .3, IX, F5.2, IX, F8. 2, IX, 

1 F8.2,1X,F8,2,1X,F5.3,1X,F5.2,1X,F8.2,1X,F8.2,1X,F8.2,1X,F5.3, IX, 

2 F5.2) 

IF ( N1-N2 ) 48,50,44 

48 IPP = IPU1 

DO 45 1= I PP, I P2 

45 WRITE (6,614) XBAS2 ( I ) , YBAS2 ( I ) , ZBAS2 ( I ) , UJ2T ( I ) , D2T ( I ) 

614 FORMAT ( 1H ,40X,F8.2, IX , F8 .2, IX ,F8 .2, 1X,F5.3, IX , F5. 2 , IX , F8.2, IX, 

1 F8.2,1X,F8.2,1X,F5.3,1X,F5.2) 

GO TO 50 
44 IPP = IPU1 

DO 46 1= I PP, I P2 

46 WRITE (6,613) XBAS1 (I) ,YBAS1 ( I), ZB ASK I ) » UJ1T ( I ),D1T( I) 

50 CONTINUE 

IF ( I HOLD 1—2 ) 90,51,51 

51 CONTINUE 

V3 = 1./VELJ3 
ZP = YJ3 
YP = -ZJ3 

WRITE (6,615) XJ3»YP,ZP»V3,DJET3 

615 FORMAT ( 1H0, 3X, 27HPR0PERT I ES OF COALESCED JET, 3X, 2HX=, F9.2, 3X, 2HY 
1,F8.2,3X,2HZ=,F8.2,3X,6HU/UJ0=,F5.2,3X,5HD/D0=,F5.2 ) 

WRITE (6,610) 

WRITE (6,616) (XBAS31 I) ,YBAS3 ( I ),ZBAS3( I ) ,UJ3( I ), D3( I), 1=1, N3) 

90 CONTINUE 

IF ( IGEOM ) 200,99,200 
200 WRITE (6,640) 

640 FORMAT (1H1) 

WRITE (6,635) 

635 FORMAT ( 1H0 , 38X ,44H*** INDUCED VELOCITIES AT CONTROL POINTS ***) 
IF (IGEOM-3) 221,221,222 

221 WRITE (6,632) 

632 FORMAT { IHO ,27X , 1HX ,8X , IHY,8X , IHZ, 12X, IHU, 14X, IHV , 14X, 1HW/ ) 

WRITE (6,634) ( XO ( I ) , YO ( I ) , ZO ( I ) ,U ( I ) , V ( I ) , W ( I ) , 1=1, NK) 

634 FORMAT ( 1H ,2 1 X ,F9. 3 , 1 X , F9.3, IX , F9.3 , 3E 15 . 5 ) 

GO TO 99 

222 WRITE (6,636) 

636 FORMAT ( 1H ,40X, 39HPRESSURE COEFFICIENTS AT CONTROL POINTS) 

WRITE (6,637) 

637 FORMAT ( 1H0,20X , 1HX,8X, 1HY,8X , IHZ* 12X, 2HCP, 14X, IHU, 14X, IHV,14X, 

1 IHW/ ) 

WRITE (6,638) ( XO ( I ) ,Y0 (I) , ZO ( I ) ,CP( I ) ,U ( I ) ♦ V ( I ) , W(I) , 1=1, NK) 

638 FORMAT ( 1H , 14X ,F9. 3 , IX , F9.3 , IX,F9.3,4E15 .5 ) 

99 CONTINUE 

RETURN 

END 
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SUBROUTINE TRANS1 (MULT, ALFA, BETA, PS ID J 


TRANSFORMS INPUT COORDINATES TO PROGRAM COORDINATES (FIXED) 
CONVERTS ANGLE OF ATTACK AND SIDESLIP TO FRSTRM DIRECTION COS. 

DIMENSION XJKIO) V YJ1(10) , ZJ1 ( 10) , DJET1 ( 10) , VELJ 1 ( 10 ) 

DIMENSION XJ2(10) , Y J2 ( 10) , ZJ2 { 10 ) , DJET2 ( 10 ) ♦ VEL J2 ( 10 ) 

C0MM0N/BLK7/ALFQ,BETQ,GETQ,Fi,F2,F3,VK0NST 
COMMON/BLKiO/XJl *YJ1,ZJ1» DJETl* VELJl,XJ2» YJ2, ZJ2, DJET2, VEL J2 

DIMENSION PSID( 1 ) 

A = ALFA*. 0174533 
B = BETA*. 0174533 
ALFQ * COS( A ) *COS( B ) 

BETQ = SIN( A) *COS(B) 

GETQ = SIN(B) 

YS = YJl(l) 

YJ1(1) = ZJ1(1) 

Z J1C 1 ) * -YS 
PSID(l) = -PSID(l) 

IF I MULT-2) 5,4,4 
YS = YJ2 ( 1 ) 

YJ2I1) * ZJ2 ( 1 ) 

ZJ2I1I * -YS 
PSIDI2I = -PSID(2) 

CONTINUE 

RETURN 

END 


SUBROUTINE TRANS2 (Y,Z,NO) 


TRANSFORMS INPUT COORDINATES TO PROGRAM COORDINATES (FIXED) 

DIMENSION Y(l),Z(l) 

DO 1 1*1, NO 
YS = Y ( I ) 

Y(I) = ZU) 

Z(I) * -YS 

RETURN 

ENO 


SUBROUTI NE TRANS3 ( Y,Z,V,W,NO,KS,NS,TNEG,NLST) 

TRANSFORMS PROGRAM COORDINATES (FIXED) TO OUTPUT COORDINATES* 
JET CENTERLINE AND CONTROL POINT COORDINATES ARE AFFECTEO 


DIMENSION 

DIMENSION 


XI ( 11,10) ,Z1(U ,10) ,UJ1 ( 11,10 ),D1( 11, 10 ) , DXDZ 1 ( 11, 10) 
X2( 11,10) ,Z2( 11,10) ,UJ2( 11,10 ),D2( 11*10 ),DXDZ2( 11,10) 



DIMENSION X1T(100),Z1T(100),UJ1T(100),D1T(100),OXDZ1T( 100) 
DIMENSION X2T < 100 ) » Z2T ( 100 ) ,U J2T ( 100 ) , D2T ( 100 ) » DXDZ2T ( 1 00 ) 
DIMENSION XBASK100) ,YBAS1 (100) ,ZBAS1( 100 ) 

DIMENSION XBAS2(100) »YBAS2(100) ,ZBAS2( 100) 

DIMENSION XBAS3U00) *YBAS 3(100 ) , ZB AS3( 100) 

DIMENSION XJ1 (10) ,YJi(10) , ZJ1 ( 10 ) , DJET1 ( 10 ) ♦ VEL J 1 ( 10 ) 

DIMENSION X J2 (10) * Y J2 ( 10 ) , ZJ2 < 10 ) , DJET2 ( 10 ) , VEL J2 ( 10 ) 

C0MM0N/BLK2/X1,ZI,UJ1,D1* DXOZ 1,X2,Z2,UJ2*D2,DXDZ2 
C0MM0N/8LK3/X1T,Z1T,UJ1T,D1T,DXDZ1T,X2T,Z2T,UJ2T,D2T,0XDZ2T 
C0MM0N/BLK5/XBAS1,YBAS1,ZBAS1,XBAS2,YBAS2,ZBAS2,XBAS3, YBAS3,ZBAS3 
C0MM0N/BLK9/MULT , IH0LD1 *KOUNT 1,I0NE,ITW0, ITHR, N 1, N2, N3* IFIXl 
C0MM0N/BLK10/XJ1,YJ1,ZJ1,DJET1,VELJ1,XJ2,YJ2*ZJ2,DJET2, VELJ2 

DIMENSION Y(l),Z(i) ,V(I),N(l) 

DO 1 1=1 » NO 
YS = Y ( I ) 

Y(I) = -Z(I) 

Z(I) = YS 
VS = VII) 

VU) = -W(I ) 

1 W( I) = vs 

KTR1 = 0 
KTR2 = 0 

IF (MULT-2) 8,5,5 

5 IF ( TNEG ) 6,6,7 

6 KTR2 = K0UNT1 
GO TO 8 

7 KTR1 = K0UNT1 

8 CONTINUE 

UJ IT ( l) = UJl(l,i) 

D1T(1) = 01(1,1) 

JK = 1 

DO 20 1=1, KS 
IF U-l) 10,10,11 

10 NF = NS+1-KTR1 
GO TO 15 

11 IF (I-KS) 12,13,13 

12 NF = NS+i 
GO TO 15 

13 NF = NLST 
15 CONTINUE 

DO 20 J=2,NF 
JK = JK+1 

UJIT(JK) = UJ1 ( J » I ) *VEL J1 ( I ) /VELJ 1 ( i ) 

20 DIT(JK) = D1(J,I)*DJET1(I)/DJET1(1) 

IF (MULT-2) 50,25,25 
25 UJ2TU) = UJ2 (1,1) 

D 2 T ( 1 ) = 02(1,1) 

JK = 1 

DO AO 1=1, KS 
IF (1-1) 30,30,31 
30 NF = NS+1-KTR2 
GO TO 35 
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31 IF U-KS) 32,33*33 

32 NF * NS+l 
60 TO 35 

33 NF = NLST 
35 CONTINUE 

00 40 J=2 ,NF 
JK = JK+1 

UJ2T(JK) = UJ2( J,I)*VELJ2(I)/VELJ2(1) 
40 D2TUK) = D2( J, I ) *DJET2 ( l ) /DJET2 < 1 ) 

50 CONTINUE 

DO 52 1=1, N1 
YS = YBASi ( I ) 

YBASKI) = -ZBASKI) 

52 ZBASKI) = YS 

IF CN2) 60,60,53 

53 00 54 1=1, N2 
YS = YBAS21I) 

YBAS2 ( I ) = -ZBAS21I) 

54 ZBAS2U) = YS 

IF <N3) 60,60,55 

55 00 56 1=1, N3 
YS = YBAS3 ( I ) 

Y8AS3U) = -ZBAS3 ( I ) 

56 ZBAS3CI) = YS 
60 CONTINUE 

RETURN 

ENO 


SUBROUTINE OUTPT < XL , ZL *DXDZ, CF , DJ ,XJ, Y J, Z J, XB* YB, ZB, VX* VY, VZ ) 

TRANSFORMS LOCAL COORDINATES TO PROGRAM COORDINATES (FIXED) 

DIMENSION CFO, 3) 

PHI = ATAN(DXDZ) 

VXT = SINIPHI) 

VYT = 0. 

VZT = COStPHI ) 

CALL ROTATE ( VX , VY, VZ ,CF, VXT, VYT, VZT, 1 ) 

CALL ROTATE (FX,FY,FZ,CF,XL»0.,ZL,1) 

XB = FX*0 J+X J 
YB = FY*OJ+YJ 
ZB = FZ*D J+Z J 
RETURN 
END 


SUBROUTINE 0UTPT1 ( XL, ZL , DXDZ,CF,KS,DJ,X J, YJ, ZJ, XB, YB, ZB, 

1 VX, VY , VZ ) 

TRANSFORMS LOCAL COORDINATES TO PROGRAM COORDINATES (FIXED), FOR 
THE SEGMENTED JETS 

DIMENSION CF ( 3 , 3 ,10 ) 
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PHI * ATAN(OXDZ) 

VXT = S I N ( PH I ) 

VYT = 0. 

VZT = COS ( PHI ) 

CALL ROTATE ( VX , VY,VZ ,CF II ,1 , KS ) » VXT* VYT, VZT* 1 ) 
CALL ROTATE { FX ,FY ,FZ ,CF ( I ♦ 1 , KS ) ,XL,0. , ZL • l ) 

XB « FX*D J+X J 
YB = F Y*0 J+Y J 
ZB = FZ*D J+Z J 
RETURN 
ENO 


SUBROUTINE PLANE (CFN1 , CFN2,CFN3,X1 , YI , Z1 , CSN 1, CSN2, CSN3, XL 1, XL2 t 
1 XL3 ,C00R1,C00R2,C00R3 ) 

COMPUTES INTERSECTION OF A GIVEN PLANE WITH A LINE 

DIMENSION CFN(3) ,CSN ( 3 ) «XL(3) * COOR ( 3 ) 

CFN< I ) = CFN1 
CFN1 2 ) = CFN2 
CFN< 3 ) = CFN3 
CSNU) = CSN1 
CSN ( 2 ) = CSN2 
C SN( 3 I = CSN3 
XL ( 1 ) = XL1 

XL < 2 ) = XL2 

XL ( 3 ) = XL3 

IL = 1 
IM = I 
IN = l 
SUBI = 0. 

IF ( ABSICSNI I) I-1.0E-04) 1,1,2 

1 IL = 0 

SUBI = CFN(1)*XL(1> 

COOR ( 1) * XL ( 1 ) 

2 IF (ABS(CSN(2) )-!• 0E-04 ) 3,3,4 

3 IM = 0 

SUBI = SUB1+CFN ( 2 ) *XL ( 2 ) 

COOR ( 2 ) = XL ( 2 ) 

4 IF (ABS(CSN(3) I-1.0E-04) 5,5,6 

5 IN = 0 

SUBI * SUB1+CFN(3)*XL(3) 

COOR I 3 ) = XL ( 3 ) 

6 0* CFN( 1 ) *X1 +CFN( 2 ) ♦Yi+CFN ( 3 )*Zl 
IF (IL+IM+IN-2) 10,30,50 

10 IF (IL) 12,11,12 

11 IF (IM) 14,13,14 

12 IP = 1 
GO TO 15 

14 IP = 2 
GO TO 15 

13 IP = 3 
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COOR IIP) * ( O-SUB i ) /CFN( I P ) 

GO TO 90 

30 IF (ID 32,31,32 

31 IP1 = 2 
IP2 = 3 
GO TO 35 

32 IF (IM) 34,33,34 

33 IP1 = 1 
IP2 = 3 
GO TO 35 

34 IP1 = I 
IP2 = 2 

35 SLOPE = CSN( IP1)/CSN(IP2) 

COOR ( I P2 ) * ( D—SUB1+CFN ( I PI )* SLOPE *XL ( IP2 )— CFN( IP l )*XL ( IP 1 ) ) / 
l (CFNI IPl )*SLOPE+CFN( IP2) ) 

COOR(IPl) * SLOPE*( COOR < I P2 )-XL (IP2) )+XL ( IP1) 

GO TO 90 

50 COEFXI = L./CSNII) 

C0EFY1 = -l./CSN(2) 

D1 = XL(1)/CSN(1)-XL(2)/CSN(2) 

C0EFX2 = l./CSNd) 

C0EFZ2 = -1./CSNI3) 

D2 = XLd)/CSN( 1 1 -XL ( 3 ) /CSN C 3 » 

CALL SOL (CFNd ) ,CFN(2) *CFN(3 ) , 0, COEFXI, COEFY 1,0. » D1 ,COEFX2,0. , 
1 C0EFZ2 ,D2«C00R ( 1 ) ,CCCR ( 2 ) , CO OR (3) ) 

90 C00R1 = COOR ( 1 ) 

C00R2 = COOR ( 2 ) 

C00R3 = COOR ( 3 ) 

RETURN 

END 


SUBROUTINE AOAMS C N, START, F INAL. H. PRINT. ICOUNT,RELB, ABSB, ISKIP, 

1 XO,XP,PAR,DDERIV) 

C 

C SUBROUTINE ADAMS SOLVES A SYSTEM OF *N* FIRST ORDER DIFFERENTIAL 
C EQUATIONS BY MEANS OF A FOURTH OROER ADAMS PREDICTOR/CORRECTOR 
C METHOD. THE STARTING SOLUTION IS BY RUNGE-KUTTA METHOD. 

C AUTOMATIC ERROR CONTROL IS OPTIONAL. 

C 

DIMENSION X< 50,5) ,VK( 50,4) »F( 50,5) ,E(50) 

DIMENSION XP(1),X0(1), PAR ( 1 ) 

C 

IBCOL = 0 

IF (PRINT) 20,10,20 
10 IF (ICOUNT) 20,31,20 
C 

20 CONTINUE 
C20 WRITE (6,400) ID,N 
I BOOL = 1 

C400 FORMAT ( 17H0PR0BLEM NUMBER I 10»5XI2HS0LUT ION OF 
C 1 13, 5X35HFIRST ORDER DIFFERENTIAL EQUATIONS.) 

C 

C SETUP INITIAL VALUES 
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00 30 I = l »N 
XII, 1) = XOU) 

30 CONTINUE 

31 CONTINUE 

IF (ICOUNT) 40,35*40 
35 ICOUNT = 9999 
40 I TEMP = 0 

BOUND * START+PRINT 
T = START 

IF { I SKI P) 45,50*45 
45 I A = 2 

IB » 4 
GO TO 2222 

50 RLTEST = 14.2*RELB 
ABTEST = 14.2 *ABSB 
FACTOR = RELB/ABSB 
BLB = RLTEST/200.0 
H = 2 • 0*H 

RUNGE-KUTTA STARTING METHOD 

1111 IA = 2 
IB = 2 
C 

2222 DO 90 J=IA,IB 

CALL DDERIV ( T, X { 1 , J-l ) ,F ( 1 , J-l ) , PAR ) 

DO 60 1=1. N 
VKU.I) = H+FU.J-l) 

X(I,J) = XII * J-l I 5*VK ( I .1) 

60 CONTINUE 

TTEMP = T+.5*H 
C 

CALL DDERIV ( TTEMP, X 1 1 , J ) ,F< 1 , J ) , PAR » 

DO 70 1=1, N 

VK 1 1 ,2 ) = H*F 1 1 , J) 

XII, J) = XII , J-l) + »5*VKI 1,2) 

70 CONTINUE 
C 

CALL DDERIV I TTEMP ,X ( 1 , J ) ,F I 1 ♦ J ) ,PAR) 

DO 80 1=1, N 
VK (1,3) = H*F 1 1 , J ) 

XII, J) = XI I, J-ll+VKI 1,3) 

80 CONTINUE 
T = T+H 
C 

CALL DDERIV ( T ,X 1 1 , J) ,F ( 1 , J ) , PAR) 

DO 85 1=1, N 

VKC I ,4 ) = H*F ( I , J ) 

XII, J) = X I I , J- 1 ) + . 16666667*1 VK ( I , 1 ) -t-2 . 0* I VK I 1 , 2 ) + 
1 VKI 1,3) )+VK(I ,4) ) 

85 CONTINUE 
90 CONTINUE 


107 



r> n o o non o o oo o oooo 


IF (IB-2) 150,3333*150 
3333 00 100 1=1, N 

XPd) = XU, 2) 

100 CONTINUE 


XP(I ) =D0UBLE INTERVAL RESULT TO BE USED IN ERROR 
ANALYSIS 

T = T-H 
H = • 5*H 

IF ( I BOOL ) 120,125,120 
120 CONTINUE 
120 WRITE ( 6,410) H 

410 FORMAT (34H0IN THE FOLLOWING CALCULATIONS H =E14.8) 

125 IF (H-. 0000001) 130,130,140 
130 WRITE (6,420) 

420 FORMAT ( 1H0, 10< 1H*) ,//// 

1 49H0EQUATI0NS CAN NOT BE SOLVED FURTHER WITHIN GIVEN 

2 14H ERROR BOUNDS.) 

RETURN 

140 IB * 3 

GO TO 2222 

150 IF (IB-3) 200,160,200 

IS ACCURACY CRITERION MET 

160 J = 3 

4444 DO 190 1=1, N 

E( I ) = ABS(XPt I J — X < I , J ) ) 

IF(E( I )-ABS( X(I,J)*RLTEST) ) 170,175,175 
170 E( I)=E( I)/ABS(X( I,J) ) 

GO TO 190 

175 IF (E( I J-ABTEST) 180,185,185 
180 E(I) = E ( I ) *FACTOR 
GO TO 190 

185 T =T-H 

IF (J-5) 3333,187,3333 

187 DO 188 K=1,N 

188 X( K , 1 ) = X(K,4) 

GO TO 1111 

190 CONTINUE 

IF (J-5) 195,6666,195 
195 IA = 4 
IB = 4 
GO TO 2222 

SHOULD ANY OF THE STARTING VALUES BE PRINTED OUT 
200 T = T-3.0*H 
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00 250 J=2*4 
T = T+H 

1 TEMP = ITEMP+1 

IF (PRINT I 210,230,210 
210 IF (T-B0UND) 230*220,220 
220 BOUND = B0UND+PRINT 
9999 CONTINUE 

C9999 WRITE (6,430) T, ( I »X( I , J) , 1=1 ,N) 

C430 FORMAT (4H0T =E14.8/ 5( 2H X, 12 , 1H=1PE12. 5 ) ) 

ITEMP = 0 
C 

230 IF ( ITEMP-ICOUNT) 240,9999,240 
240 IF (T-(FINAL-H/10.0) ) 250,999,999 
250 CONTINUE 

BEGIN ADAMS METHOD 

5555 CALL DDERIV ( T,X ( 1 ,4 ) ,F ( 1 ,4) * PAR) 

DO 260 1*1, N 

XP(I) = X ( I ,4) + .04 1666667 *H* ( 55 *0*F ( I ,4 )-59.0*F ( I, 3) 
l ♦37.0*F(I,2)-9.0*F< 1,1)) 

260 CONTINUE 

T = T+H 

CALL DDERIV ( T ,XP ( 1 ) , F ( l , 5 ) , PAR ) 

DO 270 1=1, N 

X ( 1,5) = X( I ,4)+,041666667*H*(9.0*F( I,5)+19.0*F( 1,4)- 
l 5.0*F( I,3)+F( 1,2) ) 

270 CONTINUE 

IF (ISKIP) 6666,280,6666 
280 J = 5 

GO TO 4444 

6666 IF (T-lFINAL-H/10.0) ) 295,290,290 
290 J = 5 

GO TO 999 

295 DO 300 1=1, N 

X( I ,4 ) = XU, 5) 

00 300 J=2 ,5 
F ( I , J— 1 ) = F ( I , J ) 

300 CONTINUE 

ITEMP = ITEMP+1 

TEST WHETHER COMPUTED VALUES SHOULD BE PRINTED 

IF (PRINT) 310,330,310 
310 IF ( T-( B0UND-H/10«0) ) 330, 320, 320 
320 BOUND = BOUND+PRINT 
7777 J = 4 

C WRITE (6,430) T , ( I ,X ( I , J) , 1 = 1 ,N ) 
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330 

340 

C 

C 

C 

350 

355 

C 

358 


360 

362 


365 

380 

382 


C 

999 

C999 

C440 


385 

C 


C 

C 

c 

c 

c 

c 


I TEMP = 0 

IF ( I TEMP- 1 COUNT) 340,7777,340 
IF ( I SKIP ) 5555,350,5555 

TEST WHETHER INTERVAL CAN BE DOUBLED 

DO 355 1*1, N 

IF (E(I)-BLB) 355,355,5555 
CONTINUE 

IF (PRINT) 358,380,358 
D1 = PRINT/(2.0*H) 

D 1 I =A8S ( FLOAT ( I FIX ( D1 ) )-Dl ) 

IF (Dll-.l) 362,362,360 
IF (DII-.9) 5555,362,362 
D2 = ( BOUND- T ) / ( 2 .0*H ) 

D2I=ABS ( FLOAT ( I FI X ( D2 ) ) -D2 ) 

IF TD2I-.il 380,380,365 
IF ( 021— .93 5555,380,380 
DO 382 1*1, N 
X< 1,1) = X( I ,4 ) 

CONTINUE 
H = 4.0*H 
GO TO 1111 

CONTINUE 
WRITE (6,440) 

FORMAT ( 20H0F INAL T AND XPO...) 

DO 385 1*1, N 
XP(I) = XU,J) 

CONTINUE 
FINAL = T 

WRITE (6,430) T,(I ,XI I , J) , 1*1 ,N) 

RETURN 

END 


SUBROUTINE CFCAL ( ALFQ,BETQ,GETQ,CXJ,CYJ,CZJ,CF ) 

COMPUTES DIRECTION COSINES FOR THE LOCAL COORDINATE SYSTEM, X IN 
DIRECTION OF FREESTREAM,Y NORMAL TO FREESTREAM AND INITIAL JET 
DIRECTION, Z IS XCROSSY 

DIMENSION CFO, 3) 

CF ( 1 , 1 ) * ALFQ 
CF ( 1 , 2 ) * BETQ 
CF ( 1 , 3 ) = GETQ 

CALL XPROD (CXJ,CYJ,CZJ,CF(l,l),CF(l,2),CF(l,3),CF(2,l),CF(2,2), 
1 CFO, 3) ) 

CALL XPROD (CF(l,l),CF(l,2),CF(l,3),CF(2,l),CF(2,2),CF(2,3), 

1 CF ( 3 , 1 ) ,CF(3,2),CF(3,3)) 

RETURN 

END 
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SUBROUTINE CFCAL1 < ALF, BET, GET , CXJ ,CYJ, CZ J,CF,K ) 


GIVES DIRECTION COSINES FOR THE JET-ORIENTED COORDINATE SYSTEM 
X-AXIS IS IN DIRECTION OF FREESTREAM*Y IS NORMAL TO THE PLANE 
DEFINED BY THE FREESTREAM AND JET EXHAUST DIRECTIONStZ AXIS IS 
X-CROSS-Y. CFCALl SAME AS CFCAL EXCEPT FOR PARAMETER K 

DIMENSION CFO, 3. 10) 

CF ( 1 , 1 ,K ) = ALF 
CF ( 1 ,2,K) = BET 
CF ( 1 , 3 ,K ) = GET 

CALL XPROO (CXJ,CYJ,CZJ,CF(1,1,K),CF(1,2«K),CF(1,3,K),CF(2,1,K), 
1 CF(2,2,K),CF(2,3,K) ) 

CALL XPROD (CF(1,1,K) ,CF( l,2,K) ,CF ( 1 , 3, K) , CF ( 2, 1, K ) , CF ( 2, 2, K ) , 

1 CF<2,3,K),CF(3,l,K),CF(3,2,K)yCF(3,3,K)) 

RETURN 

END 


SUBROUTINE ROTATE ( A, B, C, CF, S ,T , U, L ) 

L=0 ROTATES A,B,C INTO S,T,U, (FIXED COORDINATES TO ROTATED) 
L=1 ROTATES S,T,U INTO A, B,C, (ROTATED COORDINATES TO FIXED) 

DIMENSION CFO, 3) ,D(3) ,V(3) 

IF (L) 1,1,2 

1 D( 1) = A 
D ( 2 ) = B 
0(3) = C 
GO TO 3 

2 0(1) = S 
D ( 2 ) = T 
0(3) = U 

3 CONTINUE 
DO A 1=1,3 

A V( I) = 0. 

DO 5 1=1,3 
DO 5 J=1 ,3 
IF (L) 9,9,10 
9 M = I 
N = J 
GO TO 5 
10 M = J 
N = I 

5 V ( I ) = V( I)+0( J)*CF(M,N) 

IF (L) 6,6,7 

6 S = Vll) 

T = V ( 2 ) 

U = V ( 3 ) 

GO TO 8 
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A * V(l) 
8 = V ( 2) 
C * VO) 
8 CONTINUE 
RETURN 
END 


SUBROUTI NE XPROO < ALFI « BET1 ,GET1 » ALF2* BET2.GET2, ALF3*6ET3* GET3 ) 

COMPUTES CROSS PROOUCT OF TWO VECTORS* RETURNS A UNIT VECTOR 

ALF3 * BET1*GET2-BET2*GET1 
BET3 = ALF2*GET1-ALF1*GET2 
GET3 = ALF1*BET2-ALF2*BET1 
OENOM = SORT ( ALF3*ALF3+BET3*BET3+GET3*GET3 ) 

ALF3 = ALF3/DEN0M 
BET3 = BET3/DEN0M 
GET3 = GET3/DEN0M 
RETURN 
END 


SUBROUTI NE SOL ( All , AI2 » A13* AKI « A2 1* A22* A23* AK2, A3I* A32 » A33* AK3* 


SOLVES A SET OF THREE EQUATIONS BY METHOD OF DETERMINANTS 

DELTA = A1I*< A22*A33-A23*A32)+A21*(A32*A13-A12*A33) 

I +A31 *( A12*A23-A13*A22 ) 

XI * (AKl*(A22*A33-A23*A32)+AK2*(A32*A13-A12*A33) 

1 +AK3*( A12*A23-AI3*A22 )) /DELTA 

X2 * I All*( AK2*A33-A23*AK3)+A21*( AK3*A13-AKI*A33) 

1 +A3l*( AKl*A23-A13*AK2 ) ) /DELT A 

X3 = ( All *( A22*AK3-AK2*A32 )+A21*(A32*AKI-AI2*AK3) 

1 +A31* ( A12*AK2— AK1*A22 )) /DELTA 

RETURN 
END 
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